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ABSTRACT 

We use a global magnetohydrodynamic simulation of a geometrically thin accretion disk 
to investigate the locality and detailed structure of turbulence driven by the magnetorota- 
tional instability (MRI). The model disk has an aspect ratio H/R^ 0.07, and is computed 
using a higher-order Godunov MHD scheme with accurate fluxes. We focus the analysis 
on late times after the system has lost direct memory of its initial magnetic flux state. 
The disk enters a saturated turbulent state in which the fastest growing modes of the 
MRI are well-resolved, with a relatively high efficiency of angular momentum transport 
((a)) ^ 2.5 X 10~^. The accretion stress peaks at the disk midplane, above and below 
which exists a moderately magnetized corona with patches of superthermal field. By ana- 
lyzing the spatial and temporal correlations of the turbulent fields, we find that the spatial 
structure of the magnetic and kinetic energy is moderately well-localized (with correla- 
tion lengths along the major axis of 2.5H and 1.5H respectively), and generally consistent 
with that expected from homogenous incompressible turbulence. The density field, con- 
versely, exhibits both a longer correlation length and a long correlation time, results which 
we ascribe to the importance of spiral density waves within the flow. Consistent with prior 
results, we show that the mean local stress displays a well-defined correlation with the 
local vertical flux, and that this relation is apparently causal (in the sense of the flux stim- 
ulating the stress) during portions of a global dynamo cycle. We argue that the observed 
flux-stress relation supports dynamo models in which the structure of coronal magnetic 
fields plays a central role in determining the dynamics of thin-disk accretion. 

Keywords: accretion, accretion discs - (magnetohydrodynamics) MHD - instabilities 



1 INTRODUCTION 

Turbulence associated with the non-linear evolution of the mag- 
netorotational instability (MRI: Balbus & Hawley 1991, 1998) 
provides the dominant source of angular momentum transport 
in ionized accretion flows. This turbulence has been explored 
using numerical simulations in both the local approximation, 
in which a small patch of accretion disk is considered, and the 
global case, where a large portion of the disk is evolved. All 
of these studies have revealed that the MRI leads to sustained, 
turbulent angular momentum transport in the outward radial 
direction. 

Despite this progress, basic questions about the physics of 
the MRI and the nature of turbulence within accretion disks re- 
main unanswered. The initial numerical work on the MRI fo- 
cused on local shearing box simulations within domains that 
were a small multiple of the disk scale height H (Hawley et al. 
1995; Brandenburg et al. 1995; Stone et al. 1996) and were car- 
ried out at moderate resolutions in which the dissipation scale 
was generally much larger than in actual disks. These experi- 



ments left unresolved the issue of whether physics on the small- 
est or largest scales is important. On small scales, the extent to 
which the plasma's viscosity and resistivity (or their ratio, the 
magnetic Prandtl number) influence the large-scale dynamics of 
accretion is unclear. Resistivity - and other non-ideal MHD ef- 
fects - are indubitably important in weakly ionized protoplane- 
tary disks (Gammie 1996; Simon et al. 2011), and there is some 
evidence that a dependence on the Prandtl number persists even 
when the dissipative scales are much smaller than the disk scale 
height (Fromang & Papaloizou 2007; Lesur & Longaretti 2007; 
Simon & Hawley 2009; Longaretti & Lesur 2010). 

On larger scales, a critical question is whether global effects 
qualitatively alter the dynamics of accretion mediated by the 
MRI. It is not obvious that they should. The MRI is a local insta- 
bility, and for geometrically thin accretion disks (with H/R<. 1) 
the fastest growing linear modes are not of global scale. (We 
exclude here radiatively inefficient flows with H ^ R, for which 
the importance of global effects is self-evident.) Nonetheless, at 
least three distinct avenues via which global dynamics might af- 
fect the evolution of thin magnetized accretion disks have been 
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proposed. First, it is possible that long-wavelength fluid or mag- 
netic modes, associated with the non-linear state of the MRI it- 
self, could be important in determining either the saturated level 
of the stress, or some other property of the accretion flow. Exist- 
ing numerical investigations into this issue have yielded mixed 
results. Guan et al. (2009), using local simulations, showed that 
the primary fluid variables (B, v and p) de-correlated on a scale 
that was significantly smaller than H, supporting a fundamen- 
tally local picture of MRI-driven turbulence. Nelson & Gressel 
(2010), on the other hand, found that some secondary prop- 
erties of the turbulence (specifically, the amplitude of density 
fluctuations) converged only for domains whose size - in their 
simulations 4H x 16H x 2H - approached global scales. Second, 
additional instabilities that are not present (or only slowly grow- 
ing) on small scales could become important globally. Tout & 
Pringle (1992), for example, proposed a semi-analytic dynamo 
model in which the Parker instability played a key role. Since 
the Parker instability grows most rapidly on scales A » H, an 
implication of their model is that a global simulation of a thin 
disk ought to display quite different dynamics from a local sim- 
ulation of an otherwise identical system. Finally, even if the tur- 
bulence originating from the MRI is localized, the strength of 
that turbulence is known to depend upon the flux threading the 
disk (Hawley et al. 1995; Sano et al. 2004), which is a function 
of the boundary conditions for local simulations. Determining 
self-consistently the distribution of that flux (in the case where 
the disk as a whole is not threaded by a significant global mag- 
netic field), requires either large local (Guan & Gammie 2011) 
or global simulations. If the local flux is sufficiently strong, it 
is possible to envisage a limit in which the accretion stress is 
ultimately determined by the strength and connectivity of mag- 
netic fields in the disk corona (Tout & Pringle 1996; Uzdensky & 
Goodman 2008). There is some numerical evidence in favor of 
such a coronally-anchored dynamo (Sorathia et al. 2010). 

Prior global simulations have attempted to address some of 
the above questions (Armitage 1998; Hawley 2000, 2001; Fro- 
mang & Nelson 2006; Flock et al. 2010; Sorathia et al. 2010; 
O'Neill et al. 2010; Flock et al. 2011), but it remains compu- 
tationally demanding to run global simulations at a resolution 
high enough to allow fair comparison with local calculations. 
Many uncertainties thus remain. In this paper, we present a de- 
tailed analysis of the structure of the MHD turbulence realized 
in a global simulation of a geometrically thin disk. Compared to 
most prior work, the simulation we use has a moderately high 
resolution on the poloidal plane, but more importantly makes 
use of a higher-order Godunov scheme that is more accurate 
on test problems at Pxed spatial resolution. The same numeri- 
cal scheme is employed in the recent global simulations of Flock 
et al. (2011), but here we use initial conditions designed so that 
the saturated state is easier to resolve adequately. Our calcu- 
lations thus represent a step toward the elusive goal of higher 
accuracy, converged simulations of thin accretion disks (Haw- 
ley et al. 2011). Furthermore, we analyze the turbulence using 
methods closely related to those employed in recent local sim- 
ulations, thereby facilitating as close a comparison as possible. 
Our goal, then, is to elucidate the behavior of MRI-driven tur- 
bulence across a large range in spatial and temporal scales and 
make a connection between local and global calculations. 

The remainder of this work is organized as follows. In §2 
and §3, we give the details of the method that we use to inte- 
grate the equations of ideal Magnetohydrodynamics, along with 
information about the setup of the simulation presented here 



and how the simulation data was reduced and analyzed. In §4 
and §5 we describe the evolution of the disk, along with global 
measures and the spatial structure of the turbulent steady state. 
In §6 and §7, we examine the spatial and temporal structure of 
the turbulence and whether or not a local relationship between 
vertical flux and the magnetic accretion stress is at work in the 
simulation. Finally in §8, we summarize our results and point 
the way to future work. 



2 NUMERICAL DETAILS 
2.1 Integration Scheme 

We use the Pluto code for computational astrophysics 
(Mignone et al. 2007) to simulate the evolution of a geometri- 
cally thin accretion disk. Pluto implements a higher-order Go- 
dunov scheme to integrate general systems of conservation laws 
of the form: 

du 

--+V-F(U) = S(U) (1) 

a t 

where U is a vector of conserved variables, F(U) is a second rank 
tensor of fluxes derived from the conserved variables and S(U) 
are source terms. We utilize Pluto to integrate the equations of 
ideal MHD written in the Newtonian limit in spherical coordi- 
nates (r, 6, 4>): 

^ + V(pv)=0; -- + V-(sv) = 

ot at 

^^+V- (pw-BB + P*) =-pV<|) + G 

— + V - [(£ + P)v-(vB)B] = -pv V$ 

a t 

dB 

— + V • (vB - Bv) = 

o t 

where p is the density, s = P^/p^ is the entropy, v is the velocity, 
B is the magnetic field, E = Pjiy - 1) + p|vp + |B|V2 is the 
total energy, P* is a diagonal tensor, the non-zero components 
of which are the total pressure, P = P^ + |Bp/2, P^ is the gas 
pressure, y = 5/3 is the ratio of specific heats for an ideal gas 
equation of state, $ is the gravitational potential and G is a ten- 
sor containing geometric source terms appropriate for spherical 
polar coordinates (the precise form of which can be found in 
Mignone et al. 2007). 

Most geometrically thin accretion disks are optically thick, 
and require radiation MHD simulations for a fully consistent nu- 
merical treatment (e.g. Hirose et al. 2009). Since our goal here 
is to study the structure of disk turbulence itself - rather than 
any likely more subtle coupling between the disk's dynamics and 
thermodynamics - we adopt an approximate treatment of the 
energy equation that is designed to maintain a nearly isentropic 
evolution of the disk. We integrate both an entropy equation, 
which we use to determine the pressure in regions of smooth 
flow, and a total energy equation, which is used at shocks. Phys- 
ically, this means that we include energy dissipation and entropy 
generation at shocks (which are however negligible for the run 
presented here), whereas we discard magnetic or kinetic energy 
that is dissipated elsewhere. 

The Pluto code has been extensively tested for global 
simulations of the magnetorotational instability by Flock et al. 
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(2010). These authors report that in order to accurately repro- 
duce the Hnear growth stage of the instabiUty, a particular com- 
bination of algorithmic choices is necessary Specifically Flock 
et al. (2010) recommend the use of second order Runge-Kutta 
time-integration, second order spatial reconstruction, HLLD or 
Roe-type Riemann solvers and the upwind-constrained trans- 
port 'contact' method of Gardiner & Stone (2005) for calculation 
of the electromotive force's (EMF's) for the induction equation. 
Failure to follow these recommendations can lead to the pres- 
ence of instabilities within the evolution. For these reasons, we 
follow the algorithmic recommendations of Flock et al. (2010) 
in the simulations presented here, choosing to utilize the HLLD 
Riemann solver for reasons of robustness and computational ef- 
ficiency. 



2.2 Initial Conditions, Computational Grid & Boundary 
Conditions 

The physical inner boundary condition for thin accretion disks 
can be variously a stellar boundary layer, a stellar magneto- 
sphere (Pringle & Rees 1972), or the innermost stable orbit 
(ISCO) of a relativistic potential. Although in principle there are 
qualitative differences between these systems (for example, the 
resonance between the angular and epicyclic frequencies of a 
Newtonian potential is not present for a black hole), the basic 
properties of disk turbulence at larger radii are likely to be in- 
dependent of the inner boundary condition. For simplicity, we 
here use the Pseudo-Newtonian potential due to Paczyhsky & 
Wiita (1980): 



sound speed is specified by the power law: 



r- 1 



(3) 



This potential is a Newtonian approximation to the gravitational 
potential of a Schwarzschild black hole, whose key property is 
the existence of an ISCO at rj^co ~ ^r^, where = 2GM/c^ 
is the Schwarzschild radius. Circular orbits within this potential 
have a specific angular momentum i given by. 



'^kep 



r-l 



(4) 



Outside of rjgQQ, orbits converge toward the standard Keple- 
rian form, i^^^ = ^/r and so we expect that outside of r^gcQ, 
thin accretion disks (i.e. those with << rVL where is the 
gas sound speed) will have fluid elements on Keplerian orbits 
with i^^p = yfr even for the Pseudo-Newtonian potential given 
above. Inside of rj^^o stable circular orbits exist and the 
fluid plunges into the black hole. Provided that the inner radial 
boundary is placed deep enough within the potential, the fluid 
will become supersonic in the radial direction before reaching 
the inner boundary, and the exterior properties of the flow are 
independent of the precise details of the inner radial boundary 
(McBCinney & Gammie 2002) . For our simulation, this is accom- 
plished by placing the inner radial boundary at r = l.Sr^. 

The initial conditions for these simulations utilize a simple 
configuration, corresponding to a disk with constant H/R. The 
density is initialized to some constant, po in the midplane be- 
tween Srg < r < ISrg with a standard gaussian distribution in 
the vertical direction: 



p(Z) = poe 



(5) 



(6) 



where Z = rcosQ and H = cjQ is the disk scale height. The 



where R = rsin^, Cq is chosen to give the required value of 
H/R for the disk (here H/R = 0.07) and Rq is the disk inner 
edge. The initial pressure distribution is determined by assum- 
ing that the disk is isothermal such that P = cfp. The angular 
momentum distribution is specified by ij^^p. Finally the magnetic 
field is specified to be toroidal, = cons, within the region 
7.Srg < r < 12. Sr^ and |Z| < H. is normalized such that 
(P) = 2(Pg)/(|Bp) = 10 (where (Q) denotes a volume av- 
eraged quantity), as was used in Beckwith et al. (2008a). Note 
that the imposition of a moderately strong magnetic field on 
the hydrodynamic state disrupts any initial hydrostatic equilib- 
rium that we could hope to establish; it is for this reason that we 
choose such a simple hydrodynamic configuration initially as the 
disk rapidly relaxes to a new state at the beginning of the evo- 
lution. Finally, in order to seed the MRI, random perturbations 
are applied to the initial pressure distribution at an amplitude of 
10%. 

The simulation presented here utilized 288 x 128 x 96 zones 
in (r, 6,(1)) covering a domain l.Sr^ < r < 20r5, —5 < Z/H < 5 
and < (j) < n/2. We use a logarithmic grading for the mesh 
in the radial direction and apply strict outflow boundary con- 
ditions at both the inner and outer boundaries (i.e. no fluid or 
magnetic quantities are allowed to enter the domain). In the 
vertical direction, we locate the boundaries at ibSH so that we 
obtain stresses that are independent of the vertical domain size 
(Sorathia et al. 2010; Simon et al. 2011). The grading in the 
vertical direction is designed such that we have 32-zones per 
H in the region cover Z < \H\ i.e. 64-zones covering the re- 
gion —H < Z < H, with the remaining 64 zones covering the 
region Z > H. We apply periodic boundary conditions in the 
vertical direction; whilst these are physically unrealistic for the 
stratified global simulation that we present here, they have two 
advantages, namely that they guarantee that magnetic flux can- 
not enter or leave the domain through the vertical boundaries 
(Davis et al. 2010) and that they help to eliminate numerical 
difficulties in low density regions close to the vertical boundary 
(Reynolds & Fabian 2008). Finally, the use of a quarter- wedge 
in the (p domain, significantly reduces the computational cost of 
the calculation, whilst allowing us to capture the dominant non- 
axisymmetric modes within the flow (Hawley 2001). The grad- 
ing in this direction is uniform and we apply periodic boundary 
conditions at = 0, 7r/2. 



3 EVOLUTION DIAGNOSTICS 

3.1 Reduction and Analysis of Simulation Data 

Past experience (e.g. Beckwith et al. 2009) shows that it is es- 
sential to use high time resolution data in interpreting simu- 
lations of magnetized accretion disks. For the simulation pre- 
sented here, we performed complete three-dimensional data 
dumps twenty times per orbital period at the ISCO. We describe 
here the procedure for extracting physically relevant quantities 
from this data. The IDL routines used are available from the 
authors on request. 

The first step of our procedure is to convert the raw data 
output from the simulation (which is in binary format) into a 
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portable format (specifically, HDF5) that can be queried on any 
machine without prior knowledge of the structure of the data 
file. We regard such a procedure as essential for archival pur- 
poses. Quantities that we include in the HDF5 data set are: 
p,Pg, |Bp,V,B. HDF5 has the advantage that it can be read by 
many data analysis and visualization applications, and it is also 
possible to query individual chunks of data within a given file, 
rather than being restricted to reading the entire file, as is the 
case (for example) with raw binary data. 

The next step in the data reduction pipeline is to com- 
pute diagnostics that provide insight into the state of the mag- 
netized accretion flow described by the simulation. In past 
works, we have found it useful to consider volume-integrated 
quantities, (Q(t)), shell-integrated radial profiles, (Q(r, t)) and 
two-dimensional distributions on the r — (j) and r — 6 planes, 
(Q(^j0jf)) and {Q(r,6,t)), respectively. These are defined 
through: 



-I 



(Q(t))= Q(r,0,(l),t)r^smedrded(l) 



(Q(r,t)) = J Qir,e,(l),t)r^ sine dOdcj) 
(Q(r, (j), t)) = J Q(r, 6, (/), t) r sin OdO 



(7) 



(Q(r,0,t)) = J Ci{r,e,(t>,t)rd(i> 

For each quantity Q, we compute each of these diagnostics 
and store them as time-series, again in HDF5 format. In the 
cases where vertical integrals are performed, (Q(r, t)) and 
(Q(r, 0, t)), we compute the vertical integral over three differ- 
ent ranges. 



|Z| < 5H 
|Z| <H 
H < \Z\ < 5H. 



(8) 



We describe these as, respectively, the "disk+corona", the "disk 
body", and the "corona". Physical motivation for these distinc- 
tions is given below, although we note that what we call the 
corona is defined solely in geometric terms, rather than via a cut 
on the plasma p. We proceed similarly for the two-dimensional 
distributions of the quantity, Q on the r — (j) plane, in this case 
omitting the integral over (j) • We also consider the perturbations 
on a quantity, 5Q, which we calculate through: 



5Q(r, 0,0,0 = Q(r, 0,0,0 - 



{Q(r,e,t)} 
j rdcj) 



(9) 



with this definition, we have that {5Q(r,6,t)) =0. For the pur- 
poses of our discussion in this work, we have found it useful to 
compute these diagnostics for the same quantities as the full 
three-dimensional data dumps, p,Pg, |Bp, |V|, |B| (with vector 
quantities stored in component form) along with measures such 
as the mass-accretion rate p V", the angular momentum pi, the 
Alfven velocity v^, the sound speed c^, along with the Maxwell 
and (perturbed) Reynolds stresses. 

From these time-histories, it is easy to compute time- 
and spatially-averaged profiles of interesting quantities. We 
have found the time-averaged and volume-integral of a quan- 
tity, ((Q)), the time- and shell-averaged radial profile, ((Q(r))) 
and the time- and disk surface area-averaged vertical profile. 



((Q(0))) particularly useful, where: 



((Q(r))) ^ 



((Q(0))} = 



1 f,'jQir,t))dt 
AT jj^r^smed0d<j> 

j^jJ^J^mrMdrd^ 
AT fj-rdrd4> 



(10) 



where AT = — Tj is the time-interval over which to compute 
the time-average. 



3.2 Calculation of Power Spectra and Correlation 
Functions 

We calculate the time-averaged toroidal power spectrum, 
{{\P^(r,k,r)): 



{5Q(ir,(t),t))e-'^^'^d(t) 



(11) 



dt 



where 5Q{r,(i),t) = (Q(r, 0, t) - (Q(r, t))) / (Q(r, 0). We will 
often plot this quantity in terms of ky = Inn/ /S.y where 
Ay/H = (H/R)~^. This power spectrum can be radially- 
averaged as 



1 
AR 



PQ(r,k^))\ )dr (12) 



to 



where AR = ^2 ~ ^i- Finally we can use 
probe the ratio of power on different spatial scales by computing 



rk 



(K<^:W)) 

{5Q(ir,k^,t))e'^^Uk 



(13) 



f^' {5Q(r,k^,t))e^^^Uk^ 



dt 



where k^^^^ is the Nyquist critical frequency, /cq = 2/n is the 
largest scale in the toroidal domain, k^ is some (specified) break 
frequency and (5Q(r, /c^, t)) = j {5Q(r,(l),t)) e'^^'^'^dcl). 

The calculation of power spectra and correlation functions 
from global simulations is well-defined, but their interpretation 
is complicated by the presence of radial gradients that vary with 
time. In particular, if we simply took the Fourier Transform of 
a quantity in the radial direction without accounting for these 
gradients, then we would essentially be performing a convolu- 
tion of a flat-window function with a non-periodic function. As 
a result, the radial power spectrum would contain significant 
small scale power in order to account for radial gradients. These 
considerations would then can cause problems when it comes 
to comparing the results of global simulations against local re- 
sults. We therefore make use of a well-defined (but non-unique) 
formalism that allows us to remove (time-varying) gradients on 
large radial (and in some cases vertical) scales. 
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In developing this formalism, we have had to make sev- 
eral compromises due to the limited computational resources 
available for performing this analysis. Ideally, many of these di- 
agnostics should be computed using four-dimensional Fourier 
Transforms (three spatial dimensions + 1 time dimension) on 
individual vector components. However, since the typical dimen- 
sions of the data set to be Fourier transformed are 100^ x 1000, 
this is impractical both in terms of the memory required to store 
the data set and the time required to compute the Fourier trans- 
forms. Instead, we have had to utilize spatially-averaged scalar 
data (e.g. |Bp) in these calculations and restrict our consider- 
ation to correlation functions that either contain information 
about three-spatial dimensions, or two-spatial dimensions and 
time. 

The first set of power spectra utilized in the discussion of 
§6 are two-dimensional Fourier transforms of the magnetic and 
kinetic energy densities on the poloidal plane. These are cal- 
culated from toroidally- averaged data, (Q(r, 6,t)), weighted by 
the area element rsin^drd^. We denote the volume weighted 
quantity resulting from this procedure by (Qy(r, 0,t)). In the 
discussion of §6, we compute two-dimensional power spectra of 
this data over some restricted range in radius, < r < r2 and 
the entire vertical domain, —SH < Z < SH. Since the simula- 
tion presented here utilizes a graded mesh in (r, 0), it is nec- 
essary to first map data to a uniform mesh, which is defined 
over the range < r < —SH < Z < SH. The coordinates 
spacings, (Ar)^^^, (Az)^^^ for the uniform mesh are set so that 
(Ar),,^ < (Ar)^^„, (A^),,^ < (Az)^^„ where (Ar)^^„, (A^)^^^ 
are the smallest coordinate spacings in the simulation mesh in 
the range < r < —SH < Z < SH . We also require that 
there is an even number of cells in each dimension on the uni- 
form mesh in order to minimize the time taken by the Fourier 
transform procedure. To map the simulation data onto the uni- 
form mesh, we construct a Delaunay triangulation of the (r,6) 
mesh used in the simulation and then use this triangulation to 
perform a linear interpolation of simulation data onto the uni- 
form mesh. We denote the data mapped onto the uniform mesh 
via this procedure {Qy(r,6 ,t)) 

The next step in this procedure is to obtain a normalized 
measure of the spatial fluctuations in (Qv(r, 6, t))^ : 



^(r,0,t) = 



(14) 



This step is essential as it accounts for both radial and verti- 
cal gradients in (Qy(r, 0, t))^^^. Without this element of the 
procedure, application of a Fourier Transform in the radial di- 
rection would result in significant enhancements to small scale 
power in order to account for aperiodicity in these coordi- 
nates. Schnittman et al. (2006) describe a procedure to ac- 
complish this renormalization where the mean at each radius 
is substracted and then the fluctuations are normalized so that 
J ^^(r, 0,t)d0 = 1. This approach to renormalization of the 
fluctuations has the disadvantage that it removes power from 
all radial length scales. Our approach is to instead calculate 
((Qv(^, ^5 f))reg) t>y constructing a two-dimensional polyno- 
mial approximation to (Qv(r, 6, t))^^^ at each time t, ensuring 
that there are no radial or vertical gradients in ^(r,6,t) that 
could influence the shape of the power spectrum in either of 
these directions. 



We then take a two-dimensional Fourier transform of 



(15) 



We will often plot this quantity in terms of = 2nl/Ax 
and k^ = 2nl/Az, where Ax/H = 2(r2 - r^Xr^ + 
r2)"HH/K)"^cos(7r/4) and Az/H = 10. The time- and shell- 
averaged power spectrum ((|P(|fc|)P)) is then obtained from 



{p. 



imkim'dt 



(16) 



where \k\ = ^ /cj + P and (^(|/c|, t)) is averaged over shells of 
constant \k\. As emphasized above, this approach to investigat- 
ing the spatial power spectrum of the turbulence is well-defined, 
but non-unique. An alternative approach would be to calculate 
the power-spectrum in the radial direction using (for example) 
Bessel functions as a basis, avoiding the need to remove radial 
gradients on large spatial scales. We note however, that proce- 
dures equivalent to that described here were used extensively 
to investigate turbulence in global accretion disk simulations 
by (for example) Schnittman et al. (2006); Nelson & Gressel 
(2010) and so we use a similar approach here for consistency 
with prior works in the literature. 

A further probe of the structure of the turbulence is the tem- 
poral variability of fluctuations on spatial scales \k\. This mea- 
sure is calculated from the normalized measure of the spatial 
fluctuations, ^(r,6,t), via a three-dimensional Fourier Trans- 
form: 



(17) 



As above, we average 0^(k^,kQ,f) over shells of constant \k\ = 
y F + F to obtain ^(|/c|,/). We use this data to obtain esti- 
mates of the variability of Q on large (|fc|H/27r<l) versus small 
(_\k\H/2n) spatial scales by averaging over these 

ranges in |fc|-space. To improve statistics, we rebin the result- 
ing power spectra onto a grid in frequency space that is a factor 
of ten coarser than the original data. When doing so, we rebin 
the power in logarithmic units in order to avoid biasing that can 
occur when rebinning data that is characterized by (e.g.) red 
noise. 

The discussion of §6 and §7 makes use of auto- and 
cross-correlation functions. These measures are calculated us- 
ing Fourier transform techniques and so we proceed as above 
by mapping data on the graded simulation mesh in the region 
Ti < r < r2 and the vertical domain covering the body of the 
disk, —H < Z < H. We denote the volume-weighted data aris- 
ing from this procedure Qy(r,6,4>, t)^^g and the vertical integral 
of this data {Qy(r,(l),t))^^^. We again consider the normalized 
fluctuations in these quantities 



^(r,(/),t) = 



{Qy(r,0,^,tl,^) 



(18) 



((Q^(r,0,t)), 



Here, {Qy(r,e,(t),tX,g) and ( (Qv(r, (/>, t)) 
from one-dimensional polynomial 



, are constructed 
approximations in r to 
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[Qy(r,e,tX,g) and ( (Qv(r, t)),,^) at each 6 and t indepen- 
dently. These mean subtracted and normaUzed data are then 
Fourier transformed according to 



/// 



(19) 



The auto-correlation functions, (C^(Ar, A(^)) and 
{{C^(Ar, Acj) , At))) are then calculated as 

(C^(Ar,A0,At)) = 
j j \^(iK,k^J)\^e<^'''^^^'^^^')dk4k^df 

(20) 

To construct a time-average of (C^(Ar, A0, At)), we com- 
pute the correlation function over some period 5t = nP^j-ij^r = 
0.5(ri + r2) (where n = 2 typically). Denote this auto-correlation 
function measured at some time T during the evolution as 
(C^(Ar,A(/),At,r)). The time-average, ( (C^(Ar, Ac/), At)) ) 
is then computed as 



1 f^' 

((C^(Ar,A(/),At))) = — (C^(Ar,A(/),At,r))dr (21) 

Jti 

where AT = T2 — T-^ = NP^^^^^r > r2) is the time-interval over 
which the auto-correlation function is averaged, N = 2 typically 
and the orbital period is typically measured at r = ISrg. 

The cross-correlation function between two quantities, 
^i(r, (j), t) and ^2(^5 0? is calculated as 



(Q^^^(Ar,A(/),At)) = 



(22) 



J J J ^*(fc„fc^,/)^2(^r.^^ J)e<'^^^'^^^^'^dfc,dfc^d/ 

Here, ^l(k^,kQ,k^,t) denotes the complex conjugate of 
^liK, kQ,k^, t). The interpretation of ( (^C^^^^(Ar, A(/), At)) ) 
calculated in this fashion in that negative (positive) offsets in 
At represent fluctuations in ^^(r, (j), t) leading (trailing) fluctu- 
ations in «S2(r, 0, t). The time-average of (^C^^^^^Ar, A<p , At)"^ 
is calculated in the same way as for (C^(Ar, A(/), At)), i.e. 



((C^^^^(Ar,A<^,At))) = 

1 r^^ 

- (C^^^^(Ar,A<^,At,r))dT 



(23) 



where AT = 72 — = NP^^i^^r > r2) is the time-interval over 
which the cross-correlation function is averaged. We note that 
we will often plot both the auto- and cross-correlation functions 
in terms of Ax,Ay, A^ defined as discussed above and At in 
units of the orbital period at r = 0.5(ri + r2). 

In order to track the evolution of the auto- and cross- 
correlation functions in time, we use two approaches. The first, 
which we have found most useful when the correlation function 



is centrally concentrated (as is the case for the auto-correlation 
function) tracks the time -evolution of the maxima in the cor- 
relation function. At each At in the correlation function, we 
find the maximum amplitude of the correlation function asso- 
ciated with contributions from large (|A:|H/27r<l) versus small 
Qk\H/2n > 1) where \k\ = ^ F + F (where this distinction 
is made by calculating the correlation functions only includ- 
ing contributions from these scales) and plot the amplitude of 
this maxima as a function of time. The definition of the auto- 
correlation function means that this quantity must take its max- 
imum value at At = and so the width of the correlation func- 
tion in At provides an estimate of the lifetime of modes on large 
versus small spatial scales. The second approach, which we have 
found to be more useful when the correlation function does not 
have a well-defined shape (as is the case for the cross-correlation 
function) tracks the total amplitude of the correlation function 
as a function of time, e.g. 



(((c.,.,(At)))) = 
((c^^^^(Ar,A(/),At)))dArdA(/) 



(24) 



Because of the definition of the cross-correlation function, if 
( ( (^^1^2'^^^^) ) ) maximized at negative (positive) At, then 
fluctuations in ^^(r, 0, t) lead (trail) fluctuations in ^2(^j 0? f)- 



4 GLOBAL CHARACTERISTICS OF THE DISK 
4. 1 Evolution of the Disk 

The evolution of the disk proceeds in a fashion consistent with 
prior toroidal field models, e.g. Hawley (2000); Hawley & BCrolik 
(2002); Beckwith et al. (2008a). Modes with high poloidal and 
high toroidal wave numbers grow first (as described in Hawley 
et al. 1995, see top left panel of Figure 1). These modes gradu- 
ally assemble into structures characterized by smaller toroidal 
and poloidal wave numbers, as can be seen in the top right 
panel of Figure 1. This process continues for a further 2.5 or- 
bits (measured at r = 15r5, until after about three orbits at 
this radius the amplitude of the turbulent fluctuations have be- 
come sufficient to drive accretion into the central object (bot- 
tom left panel of Figure 1). By five orbits at r = 15r5, tur- 
bulent fluctuations in the magnetic field have expanded to fill 
the entire radial domain of the simulation and the disk has 
entered the quasi-stationary state (bottom right panel of Fig- 
ure 1). The simulation is evolved for a total of 20 orbits at 
r = 15r5; the data of Figure 2 shows the state of the disk at this 
time. Note the large amplitude non-axisymmetric surface den- 
sity fluctuations and equipartition-strength magnetic fields evi- 
dent in these plots. Beckwith et al. (2008a) found that models 
initialized with a net toroidal field possessed corona with mag- 
netic fields with volume-averaged strengths a factor of approxi- 
mately three below equipartition. The data presented in Figure 
2 suggest that whilst that conclusion is approximately correct 
in a volume-averaged sense; it is possible to form equipartition- 
strength magnetic fields in non-axisymmetric structures, remi- 
niscent of those suggested by Spruit & Uzdensky (2005). 

Fromang & Nelson (2006) demonstrate that net toroidal 
flux present in the initial state of their models is rapidly expelled 
from the disk into the coronal region, such that sustained MHD 
turbulence within the disk is dependent on the small scale dy- 
namo exhibited in stratified shearing boxes (see e.g. Davis et al. 
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2010; Gressel 2010; Simon et al. 2011). Figure 3 plots the time 
history of the toroidal flux, F^(t) in the simulation presented 
here, where: 



J r—5rs -J 6 



B^(ir,e,(l) = n/4,t)rdrde 



(25) 



The evolution of this quantity is similar to that reported by Fro- 
mang & Nelson (2006); the initial net toroidal flux distribution 
is expelled from the region of integration through the radial 
boundaries over the first five orbits in the outer disk (i.e. soon af- 
ter the linear growth phase of the MRI) and is replaced with flux 
of opposite sign. After approximately eleven orbits at r = ISr^, 
a small scale dynamo reminiscent of that observed in vertically 
stratified shearing box simulations produces toroidal flux of op- 
posite signs with a cycle of approximately five orbits, as pre- 
viously reported by O'Neill et al. (2010). The relatively rapid 
reversals of the toroidal flux suggest that the late-time analysis 
of the disk structure, presented below, ought to be largely inde- 
pendent of our choice of a net (rather than a zero) flux initial 
condition. Figure 4 shows the structure of ((^c/)))^ averaged 
over each period of the dynamo cycle shown in Figure 3, specif- 
ically orbits 10 - 12, 12 - 14.5 14.5 - 16 and 16 - 19 where 
the orbital period is measured at r = 15r5. These periods are 
chosen to correspond to periodic minima and maxima in the dy- 
namo cycle shown in Figure 3. The toroidal magnetic field is 
organized both radially and vertically over large spatial scales 
and the sense of this organization varies over the course of the 




Figure 2. Volume renderings of gas density (top panel) and gas (3 pa- 
rameter (bottom panel) at t = 7000. 



dynamo cycle in a quasi-periodic fashion, reminiscent of the cy- 
cles seen in vertically stratified shearing box simulations (Davis 
et al. 2010; Gressel 2010; Simon et al. 2011) and in previous 
global simulations at individual radii (O'Neill et al. 2010). 



4.2 Angular Momentum Transport and Resolved 
Turbulence 

A useful quantitative measure of the turbulence in these simu- 
lations is (a{;-''(t)) = (w^^'^Ct)) / (P,,„(t)) where w/^ = 
p5Vj.5v^, = —Bj.B^ are the Reynolds and Maxwell accre- 
tion stress respectively, W^^ = W^^ + is the total stress and 
Pg j^ are the gas and magnetic pressures. The time-evolution of 
^^/,m,f ^ is shown in the left-hand panel of Figure 5, whilst that 
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t/P„,b(r=15rs) t/P,,b(r=15rs) 



Figure 3. Time history of the toroidal magnetic ^lux, F^^, computed from 
the toroidal magnetic field evaluated at (/) = n/4 and integrated over 
the radial region 5 < r < 15 and the entire vertical domain of the simu- 
lation, |Z| <5H. 





d 0.3 




Figure 5. Time history of accretion stress normalized to the gas pres- 



, t,f,m 

sure, {oCg 



iw 



t,f,m\ 



I \Pg) (top panel) or the magnetic pressure, 

{cim'^"^ = ^^r/'^) / {Pm) (bottom panel). All quantities were inte- 
grated inside the 'disk-bod/ (|Z| < H, black lines) or the 'disk+corona' 
(red lines) over the radial region 5 < r < 15. In both panels, solid 
lines denote the total accretion stress, (^W^^^, dashed lines the Maxwell 

stress, ^^/^^ and dot-dash lines the Reynolds stress, (w^^^. Vertical 
dotted lines indicate the averaging interval. 



«B.(r,9)» X 10"^ 



Figure 4. Evolution of the azimuthal average of the toroidal magnetic 
field, ^ (^B(j) (r, 0)) ) on the poloidal plane, averaged between orbits 10 — 
12 (top left), 12 - 14.5 (top right), 14.5 - 16 (bottom left) and 16 - 19 
(bottom right) where the orbital period is measured at r = 15r5. 



of (ct{;^'^^ is shown in the right-hand panel of this figure. The 
data in these figures are computed by volume-averaging sim- 
ulation data over the radial range 5 < r/r^ < 15 such that 
we exclude regions of the disk likely to be influenced by the ra- 
dial boundaries and either the entire vertical extent of the sim- 
ulation (which we refer to as "disk+corona", denoted by red- 
lines in these figures) or ±H from the midplane (which we re- 
fer to as the "disk body", denoted by black-lines in these fig- 
ures) . The data of this figure are consistent with those of verti- 



cally stratified simulations computed at similar resolutions both 
with and without net flux configurations (see e.g. Guan et al. 
2009; Simon et al. 2011; Davis et al. 2010); specifically, we find 

((^0) ^'^■^^ ii^m)) ^^■'^ measured in the disk- 

body between t = 11.5 — 20 T^^-i^ir = ISr^) (the choice of av- 
eraging interval is such that it coincides with two dynamo cy- 
cles as detailed in the previous section). Note though that we 
find higher Reynolds stresses than reported by these authors, 
with ((«{)) = 0.3 - 0.5 ^^a^^^ compared to ((«{)) = 

0.25 ^^a^^^. The magnitude of ((^^^^^ is also significantly 
higher than reported in previous global calculations of magne- 
tized thin accretion disks. Fromang & Nelson (2006) reports 
~ 5.0 X 10~^ for a toroidal field model computed at approx- 
imately twice the resolution of the calculation here, whilst So- 
rathia et al. (2010) report for a model begun with a weak zero 
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Figure 6. Spatial dependence of number of zones per characteristic 
wavelength of the MRI, ((Qc/))) on the r — 6 plane. Simulation data 
was averaged in the toroidal direction and time-averaged over the pe- 
riod AT = 10 - 20Porb(r = 15rs). 



net poloidal flux configuration 6.0 x 10 ^ < < 9.0 x 10 ^ for 
a simulation of comparable resolution to that presented here^ 

One possible explanation of the substantially different lev- 
els of angular momentum transport measured in this work is 
that the simulation presented here has a significantly different 
resolution (measured in terms of the MRI) compared to those 
of (e.g.) Fromang & Nelson (2006); Sorathia et al. (2010). Fro- 
mang & Nelson (2006) find that for the MRI to be resolved in a 
global simulation begun with a net toroidal magnetic field, re- 
quires the fastest growing unstable mode of the toroidal field 
MRI to be resolved by at least 5 zones. These authors demon- 
strate that if this criterion is not satisfied over a significant frac- 
tion of the disk volume in a time-averaged sense, then the MRI 
becomes under-resolved and angular momentum transport de- 
cays. This requirement can alternatively be expressed in terms of 
the toroidal "quality factor" described by Hawley et al. (2011): 



(26) 



where A^^^ is the "characteristic" wavelength of the MRI 
(closely related, but not precisely equal to the wavelength of the 
fastest unstable mode, see Hawley et al. 2011). Figure 6 shows 
((Qc/,)) calculated over the interval AT = 11.5 - 20 P^rbC^ = 
ISr^). The data of this Figure demonstrate that the MRI is well 
resolved, i.e. Q^^ > 5 throughout the magnetized region of the 
disk, at least in terms of the criteria specified by Fromang & 
Nelson (2006). We also note that > 10 within the coronal 
region. The MRI is well-resolved here because the density de- 
creases due to vertical stratification. Magnetic fields rise buoy- 
antly from the disk-body into the corona at constant (if not in- 
creasing) field strength (see e.g. Figure 4). As a result, the Alfven 
velocity in the corona is (significantly) greater than in the disk 
body and hence A^^^ increases. Beyond this, buoyant magnetic 
structures tend to become dominated by power at large scales 
(e.g. Suzuki & Inutsuka 2009; Blackman & Pessah 2009). Ver- 
tically averaging inside the disk body, we find ~ 10, cor- 




^ Note that Hawley & Krolik (2002) reports 
H/R~0.15. 



0.1 for a disk with 



100 150 200 250 300 

t/Porb(risco) 



Figure 7. Time histories of the mass-accretion rate, M^rjsco^ (top 
panel) and reduced specific angular momentum flux, jneti^isco^ (bot- 
tom panel) evaluated at the ISCO. Both these quantities were evaluated 
from simulation data shell-integrated over the 'disk-body', |Z| < H. In 
both panels, solid black lines denote simulation data, vertical dotted 
lines the averaging period, red dashed lines the mean value calculated 
over this period and in the right panel, black dot-dash lines the angular 
momentum of a circular orbit at the ISCO. 



responding to a characteristic mode with toroidal wave num- 
ber m = 40. These results demonstrate that the simulation pre- 
sented here is of approximately the same effective resolution as 
the well-resolved models described in Fromang & Nelson (2006) 
which exhibited sustained accretion stresses over many hun- 
dreds of inner disk orbits. That is, the different levels of angular 
momentum transport measured in this work are not the result 
of substantial differences in the effective resolution of the simu- 
lations. We will return to the origin of this discrepancy in §5.3. 

That > 5 is a necessary, but not sufficient demonstra- 
tion that the MRI is well-resolved in the simulation presented 
here. The data of Figure 5 shows an apparent secular decrease 
of ^a^(t)^, which is not observed in ^aJJJ(t)^. Inspection of 
simulation data reveals that both that the volume integrated gas 
and magnetic pressure both decline after t ~ SP^rbCr = ISr^), 
i.e. after the linear growth phase of the MRI has completed. The 
decline of these quantities does not occur in "lock-step" how- 
ever, magnetic pressures decrease more rapidly than gas pres- 
sure and p = Pg/Pm decreases from ^ = 10 initially to ^ = 25 at 
the end of the simulation, behavior consistent with that reported 
by Beckwith et al. (2008a). The data of Figure 5 show that as 
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Figure 8. Radial profiles of characteristic velocities within the accre- 
tion flow, {{|V|(r)>>. All quantities are time- averaged over the period 
AT = 10 — 20Porb(r = ISrg) and shell-averaged over the 'disk-body' , 
\Z\ < H. Solid lines denote the accretion velocity, {{vacci^))) ^ dashed 
lines the orbital velocity, {{^4)i^)))> dot-dash lines the sound speed, 
{{cs(r))) and en-dash lines the Alfven speed, ( (v^(r)} } . The red dashed 
lines show the orbital velocity for a Keplerian angular momentum profile 
(included for reference) . Note that the flow becomes both super- Alfvenic 
and super-sonic before reaching the inner radial boundary. 



this process occurs, (^a^(t)j remains constant, suggesting that 
the MRI remains well-resolved throughout this process (see e.g. 
Hawley et al. 2011). We therefore conclude that the secular de- 
crease in ^a^(t)^ is due to evolution of the accretion flow itself, 
rather than the MRI being under-resolved here. 

A further test as to whether the MRI remains well re- 
solved is to examine fluxes of mass and angular momentum 
through the ISCO. Noble et al. (2010) demonstrated that if the 
MRI becomes under-resolved during the evolution, then the ra- 
tio of these two quantities shows a secular increase. Figure 7 
shows M(rjsco,t) — ~ {p^X^isco^^)) ^i^d the net specific an- 
gular momentum carried through this surface, jneti^isco^O — 
Urisco^ t)/M(rjsco^ 0, where 



(27) 



in analogy to the discussion of McKinney & Gammie (2002). 
These quantities are integrated within the "disk-body" (as de- 
fined above); extending the integral outside of this region in- 
cludes contributions from low density, low angular momentum 
and rapidly inflowing fluid present in the corona. The purpose 
of this diagnostic is to demonstrate that turbulent fluctuations 
within the disk body remain well-resolved and so we exclude 
the coronal material here. The data of Figure 7 demonstrate 
that over the duration of the simulation, there are no long term 
trends in j^^^ evaluated at the ISCO; we note that this quantity 
time-averaged over AT = 11. 5 — 20 P^rh^r = ISr^) is reduced 
by approximately 0.5% compared to the angular momentum of 
a Keplerian orbit at the ISCO, a result consistent with the work 
of Beckwith et al. (2008a,b) who found significant reductions 
in electromagnetic stresses at and inside the ISCO for models 
begun with toroidal as compared to poloidal magnetic field dis- 
tributions and hence net angular momentum fluxes consistent 
with that of a Keplerian orbit at the ISCO. The mass accretion 
rate through this surface displays more complex behavior with at 
least two different states evident; a high state prior to 150 ISCO 
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Figure 9. Vertical profile of turbulent velocity fluctuations, 
^ ^|V^|(z/H)^ ^, in units of the gas sound speed evaluated in the 
midplane, {{c^izlH = 0))) . Both quantities are time-averaged over the 
period AT = 10 — 20Porb(r = ISrg) and radially- averaged over the 
region S <r /rg <lb. Dashed lines denote \{z I H))) , dot-dash lines 
((|V^|(z/H))) and en-dash lines ( (|V^ |(z/H)) ). 



orbits (corresponding to 11.5 orbits at r = ISrg) and a subse- 
quent low state. Fluctuations about the mean in the low state 
occur within ~ ±1 standard deviation; suggesting that disk has 
entered a quasi-stationary state during this period. This is again 
consistent with the results of Beckwith et al. (2008a) who found 
that it took approximately 12 outer disk orbits for toroidal field 
models to establish a quasi-stationary state in terms of the mass 
accretion rate. 



5 STRUCTURE OF THE DISK 

Having demonstrated that the development of non-linear mag- 
netohydrodynamic turbulence within the disk conforms to the 
expectations of previous studies and that the MRI remains well- 
resolved throughout the simulation, we now characterize the 
radial and vertical structure of the disk. This is accomplished 
through the shell- and time-averaged radial profiles, where the 
limits of integration are restricted to the "disk-body" (i.e. |Z| < 
H) and disk area and time-averaged vertical profiles averaged 
over the radial range 5 < r/r^ < 15. Both of these diagnostics 
are calculated as described in §3.1 and time-averaged over the 
period AT = 11.5 - 20 P^,^{_r = ISrg). 
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Figure 10. Radial profiles of components of the magnetic 
field, time- averaged and vertically integrated within the disk- 
body, ^^|B^p(r/r5)^^, normalized to magnetic field strength, 
((|B|2(r/rs))). Solid lines denote ( (|Bn^(r/rs)) ), dashed lines 
((|B^|2(r/rs))) and dot-dash lines {{\B'\\r/rs))). 

5.1 Velocities 

We begin by considering radial profiles of several characteristic 
velocities within the accretion disk body, including the Alfven 
speed ((v^(r))), sound speed ((c^r))), accretion velocity, 

«v.„(r))) = - <(pv,(r))) / ((p(r))) (28) 

and orbital velocity 

((v^(r))) = r-' ((p£(r))) / ((p(r))) (29) 

Radial profiles of these quantities are shown in Figure 8. The 
data of this figure show that, well outside the ISCO (r > Sr^ 
which we term the "outer" disk) the disk is characterized by 
highly supersonic orbital motion, subsonic Alfven velocities (cor- 
responding to p = Pg/Pj^ ^ 25) and slow inward radial drift 

characterized by {{VacXr))) / {{^(t)(^))) S 10~^. As we move 
inwards through the disk towards the ISCO, the accretion ve- 
locity increases, eventually exceeding the Alfven speed at ap- 
proximately the radius of the ISCO and the sound speed just 
inside this point, i.e. the flow outside of the ISCO is upstream 
of both a sonic and Alfven point, such that the inner boundary 
should not influence the structure of the outer accretion disk 
(as described by McKinney & Gammie 2002, see §2). We also 
note that both the ordering of velocities within the outer disk 
and their radial dependence is consistent with those reported 
by Hawley (2000), with two caveats. Firstly, the accretion ve- 
locity in the outer disk is approximately an order of magnitude 
less and furthermore, the Alfven point is closer to the ISCO by 
AR ~ 1 than reported by Hawley. The simulations reported by 
Hawley are, however, moderately thick accretion torii. A bet- 
ter point of comparison then is the data of Fromang & Nelson 
(2006); Reynolds & Fabian (2008) who examine the accretion 
velocity within thin magnetized accretion disks. Fromang & Nel- 
son (2006) find radial velocities consistent with those demanded 
by the equations of standard thin accretion disk theory (Pringle 
1981), where the inflow timescale is much longer than the or- 
bital timescale. Clearly, our results are consistent with this re- 
quirement. Reynolds & Fabian (2008) demonstrate that the form 
of the accretion velocity within the disk is consistent with slow 
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Figure 11. Toroidal power spectrum of magnetic field strength, 
^ ^ I (^iBpC^y)) I ^ plotted as ^ ^ I (P|B|2(fcy)) | ^ ^ such that the 

J -axis is proportional to the fractional contribution to the total power 
per logarithmic interval in ky. In calculating this quantity, simulation 
data was first integrated vertically inside the 'disk-body' , |Z| < H, 
to yield (|B|2(r/rs, c/), t)). This quantity was then Fourier transformed 

along the (^-axis to obtain | (-PiBpC^y)) | • Finally, | (-P|Bp(^y)) | 
time- averaged over the period AT = 11.5 — 19 Pgr bi^ = l^rg) and 
radially-averaged over the region 5 < r/rg < 15. 

inward radial drift well outside of the ISCO transitioning to bal- 
listic infall close to the ISCO. Additionally, Reynolds & Fabian 
(2008) find that the Alfven point is located a distance ~ H in- 
side the ISCO, again consistent with the results presented here. 

Figure 9 shows the vertical structure of turbulent velocity 
fluctuations, {{\5V^\(z/H))) in units of the gas sound speed 
evaluated at the midplane, {{cXz/H = 0))). In calculating the 
data shown in this figure, we averaged simulation data over the 
disk surface area between 5 < r/r^ < 15 over all (j) and time- 
averaged over the period AT = 10 — 20Porb(^ = l^r^). The data 
of this figure are directly comparable with that of Figure 1 1 of 
Fromang & Nelson (2006). We find a similar structure to the 
turbulent velocity fluctuations within the disk to that described 
by these authors; fluctuations are dominated by the radial com- 
ponent of the velocity, which are approximately a factor of two 
large at the midplane than fluctuations in the (p— or 6— com- 
ponents of the velocity. Fromang & Nelson (2006) associate this 
enhanced fluctuations in the turbulent radial velocities with ra- 
dially propagating spiral density waves, which are evident in the 
simulation presented here (as can be seen from Figure 2). At 
the midplane, the radial velocity fluctuations are approximately 
10% that of the sound speed here, rising to approximately 20% 
of the midplane sound speed at |Z| = 3H, a profile closely fol- 
lowed by turbulent fluctuations in the (j) — and 6 — components 
of the velocity, but lower in amplitude at all heights by factors 
-2. 



5.2 Magnetic Fields 

Figure 10 shows the radial profile of the shell- and time- 
averaged magnetic field, ((|B'p(r)))5 (i = r, 0,0), normal- 
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Figure 12. Ratio of power in the magnetic field strengthon small toroidal 
scales (kyH/2n > 0.3, corresponding to ~ 3H) to power in this same 
quantity on large toroidal scales ikyH/2n < 0.3). In calculating this 
quantity, simulation data was integrated vertically inside the 'disk-body' 
, |Z| < H, to yield (^\B\^(r/rs,4>,t)y This quantity was filtered in 
Fourier space to include contributions from either kyH/2n < 0.2 or 
kyH/2n < 0.2 and then transformed back into real space. Finally, the 
resulting filtered data was averaged in the toroidal direction and time- 
averaged over the period AT = 11.5 — 19Porb(^ = l^rs). 



ized to the sum of these quantities at each radius. The mag- 
netic field is predominately toroidal outside of the ISCO; ap- 
proximately 10% of the total energy in the field is radial while 
vertical field accounts for only ~ 2% of the total field energy. 
Inside the ISCO, the balance in field components changes as the 
fluid plunges into the black hole; here the relative importance 
of the radial field increases with decreasing radius and there is 
a corresponding decrease in the relative importance of toroidal 
and vertical magnetic fields. This result, taken in combination 
with the data of the Figure 8 suggests that the behavior of the 
magnetic field inside of the ISCO is determined by flux-freezing, 
rather than turbulence. To confirm this suggestion, consider the 
one-dimensional toroidal power spectrum of the magnetic field 



strength, (^P\^\2(ky)j j j, shown in Figure 11, calculated 

as described in §3.2 where the radial average was computed 
over 5 < r/vg < 15. Note the location of the break in the power 
spectrum at kyH/2n = 0.3, corresponding to a spatial scale of 
3H (the significance of which will be discussed in §6). Figure 
12 shows the radial profile of the ratio of power in the mag- 
netic field strength on toroidal scales smaller than this break 
{kyH/2n > 0.3) to power in this same quantity on toroidal 
scales larger than this break (_kyH/2n < 0.3) (calculated as de- 
fined in eqn. 13 with kiH/2n = 0.3, see §3.2). In the region 
4 < r/r^ < 8, we find that the ratio of these two measures is 
~ 0.13, while inside of 4rg, power on toroidal scales smaller 
than the break drops rapidly compared to that on toroidal scales 
larger than the break. That is, as one approaches the ISCO from 
larger radii, we measure decreasing power in small scale turbu- 
lent fluctuations of the magnetic field in the azimuthal direction, 
with a corresponding increasing in power in large scale modes, 
exactly as we would expect if the dynamics of the magnetic field 
were controlled by flux-freezing, rather than small scale turbu- 
lent fluctuations. 

Figure 13 shows the vertical structure of the magnetic field 
strength, {{\B\^(z/H))) and the magnetic field components. 




Figure 13. Vertical profiles of magnetic field strength, 
((|B|2(z/H))) and ( (|B^ |2(z/H)) ) (left panel) and (3 = 
2((Pg(z/H)))/((|B|2(z/H))) (right panel). All quantities are 
time- averaged over the period AT = 11.5 — 19 Pgr bi^ = l^rg) and 
radially-averaged over the region 5<r/r5<15.In the left-hand panel, 
solid lines denote ( (|B|2(z/H)) ), dashed lines denote ( (|Bn^(z/H)) ), 
dot-dash lines ( (iB'^ |2(z/H)) ) and en-dash lines ((|B^|2(z/H))). 



{{\B^\^(z/H))), along with that of the gas p parameter, where 

l5 = 2{{p^(,z/H)))/{{\B\Hz/H))) 

As with the vertical profile of turbulent velocity fluctuations dis- 
cussed previously, to compute the data shown in this figure, we 
averaged simulation data over the disk surface area between 
5 < r/r^ < 15 over all (p and time-averaged over the period 
AT = 10 - 20Po,b(r = 15r5). The vertical profile of each of 
the magnetic field components, along with that of the magnetic 
field strength, show a similar vertical structure; approximately 
constant inside of |Z| < 2H and the falling by a factor ~ 2 — 3 
at higher latitudes. The relative strengths of the magnetic field 
components is as expected from the radial profile shown in Fig- 
ure 10, the field is predominantly toroidal, with radial contri- 
butions at roughly the 10% level and small contributions from 
vertical fields. Unlike in Fromang & Nelson (2006), there is no 
hint that the magnetic field topology changes as one moves from 
the 'disk-body' to higher latitudes, perhaps due to our utiliza- 
tion of periodic vertical boundary conditions. The vertical pro- 
file of the gas P parameter is shown in the right-hand panel 
of Figure 13. We find that this quantity varies by roughly an 
order of magnitude over the vertical extent of the disk, from 
^ ~ 30 in the disk midplane, to ^ = 3 at |Z| = 3H. On average, 
therefore, the corona is only moderately magnetized in this sim- 
ulation. Recall, however, that there are patches of the corona 
where p < 1, as evident in Figure 2, emphasizing the impor- 
tance of non-axisymmetric structures. The vertical dependence 
of p the disk magnetization intermediate between the results 
of Fromang & Nelson (2006) and Beckwith et al. (2008a). The 
former set of authors report stronger magnetizations in the disk 
midplane (p ~ 10) and slightly stronger magnetizations in the 
corona (^ ~ 1), whilst the latter set of authors report magneti- 
zations weaker in the midplane (p ~ 100) and slightly stronger 
in the corona (p ~ 1). The range in magnetization found here 
is therefore similar to that found in Fromang & Nelson (2006), 
but smaller than that of Beckwith et al. (2008a), this latter con- 
trast perhaps due to the use of full GRMHD by Beckwith et al. 
(2008a). We also note here that even though the magnetic field 
strengths found in this simulation are approximately a factor of 
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Figure 14. Radial profile of accretion stress, ((^rV^^-^)) ^^^^^^ ^^^^^ 
compared to prediction of a model that assumes zero ISCO stress 
(dashed line). Simulation data were time-averaged over the period 
AT = 11.5 — 19Porbi^ = ISrg) and shell-averaged over the 'disk-bod/ , 
\Z\<H. 



three weaker than those reported by Fromang & Nelson (2006), 
the volume averaged Maxwell stress levels found here are nearly 
an order of magnitude greater than those reported by these au- 
thors. 



5.3 Accretion Stresses 

Our next probe of the radial structure of the disk is the accretion 
stress. As we noted in the previous section, the magnetic field 
strength in the simulation presented here is a factor of three 
weaker than that reported by Fromang & Nelson (2006), whilst 
the volume averaged accretion stress is approximately an order 
of magnitude greater than reported by these authors. Taken to- 
gether, this result suggests that at fixed magnetic field strength, 
the simulation reported here exhibits accretion stresses a factor 
thirty larger than those reported by Fromang & Nelson (2006). 
Since we have already found that the effective resolution of the 
simulation presented here is the same as that found by Fromang 
& Nelson (2006) and that the MRI is well-resolved, the origin 
of this discrepancy must be related to the physical properties of 
the turbulence in this simulation. We will therefore analyze the 
properties of the accretion stress directly in order to understand 
the discrepancy. 

The data of Figure 14 shows the time-averaged radial pro- 
file of the total accretion stress, ((^rV^^^)) compared to the 
expectation from a combination of arguments regarding angu- 
lar momentum conservation and the assumption that the ISCO 
is stress free (see e.g. Hawley & Krolik 2002; Beckwith et al. 
2008b, for a detailed description) : 



m)M^(^_l,sco 



n/2 



m 



(30) 



Here, Mjg^Q is the time-averaged mass accretion rate through 
the ISCO, r2(r) is the time-averaged orbital angular velocity, 
iisco is the time-averaged specific angular momentum of the 
fluid at the ISCO and = R^n(R) is the time-averaged spe- 
cific angular momentum of the fluid. The data of Figure 14 
makes clear that angular momentum conservation combined 
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Figure 15. Toroidal power spectrum of the Maxwell stress, 



, calculated as described for the magnetic field 



strength in Figure 11. 




Figure 16. Ratio of power in the total accretion stress, ((^rV^^^)) 
on small toroidal scales ikyH/2n < 0.3, corresponding to ~ 3H) to 
power in this same quantity on large toroidal scales (kyH/2n < 0.3), 
calculated as described for the magnetic field strength in Figure 12. 



with a stress-free boundary condition applied at the ISCO pro- 
vides an excellent description of the turbulent stresses within the 
disk body for r > 4rg . As one moves in towards the ISCO from 
r = 4r5, the radial form of the total stress deviates from the pre- 
diction of the stress-free ISCO assumption and well inside of the 
ISCO, the stress profile becomes singular as the inner bound- 
ary as approached from above. The preceding discussion sug- 
gests that departures from the stress-free inner boundary con- 
dition originate from the increasing importance of flux freezing 
effects as the ISCO is approached from above. We demonstrate 
this using the same approach outlined above for the magnetic 
field strength. Figure 15 shows the toroidal power spectrum 



of the Maxwell stress. 



calculated as de- 



scribed in §3.2 where the radial average was computed over 
5 < r/r^ < 15. As was the case for the toroidal power spec- 
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Figure 17. Vertical profile of accretion stress, (^(Wj-(^(z/H)^^ , in units 
of the gas pressure evaluated in the midplane, ^ ^Pg(z/H = 0)^ ^. Both 
quantities are time-averaged over the period AT = 11.5 — 19Porb(.^ = 
ISrg) and radially- averaged over the region 5 < r/rg < 15. Solid 
lines denote the total stress, ^ ^W5^(z/H)^ ^, dashed lines denote the 

Maxwell stress, ^ (w^J^^z/H)^ ^ and dot-dash lines the Reynolds stress. 



trum of the magnetic field strength, a break in this power spec- 
trum is evident on scales kyH/2n = 0.3. Figure 16 shows the 
radial profile of the ratio of power in the total accretion stress 
on toroidal scales smaller than this break (kyH/2n > 0.3) to 
power in this same quantity on toroidal scales larger than this 
break (kyH/2n < 0.3) (calculated as defined in eqn. 13 with 
kiH/2n = 0.3, see §3.2). The data of these panels clearly illus- 
trate that contributions from small scale fluctuations to the ac- 
cretion stress decrease dramatically in importance inside of 4rs, 
confirming the suggestion that inside of 4rg, accretion stresses 
are controlled by flux freezing, rather than turbulence. 

Next, we consider the vertical structure of the accretion 
stress, ^ (w^'^'^(z/H)^ ^ normalized by the gas pressure evalu- 
ated at the midplane, (^(^Pg(z/H = 0)^Y These data are shown 
in Figure 17 and were obtained from averaged simulation data 
over the disk surface area between S < r /r^ < IS over all (j) and 
time-averaged over the period AT = 10 — 20Porb(^ = l^r^). 
Notably, the vertical profiles of the total. Maxwell and Reynolds 
stresses are all approximately constant for |Z| < 2H, a result 
consistent with those of Simon et al. (2011) for vertically strat- 
ified shearing box simulations performed at resolutions of 32 
zones per scale height in the "high" state. Note that both Fro- 



mang & Nelson (2006) and Sorathia et al. (2010) found a dis- 
tinct 'double-peaked' form to the vertical accretion stress profile, 
reminiscent of the 'low' state reported by Simon et al. (2011), 
where resistive effects had resulted in a disk that had little-to-no 
angular momentum transport and low levels of magnetohydro- 
dynamic turbulence. Intermittent "double peaked" structures in 
the vertical accretion stress profile are also evident in the space- 
time diagram of the Maxwell stress presented in Davis et al. 
(2010) Such double peaks are also evident in the data of Simon 
et al. (2011) for the "high" state. However, when these data 
are time-averaged over many local orbital periods, these dou- 
ble peaked structures are removed in the "high" state, but not 
for the "low" state. The data of both Fromang & Nelson (2006); 
Sorathia et al. (2010) was time-averaged over many orbital pe- 
riods to obtain the vertical stress profile, suggesting that both of 
these simulations were in a state equivalent to the "low" state 
described by Simon et al. (2011). 

A possible explanation as to why the simulations of Fro- 
mang & Nelson (2006); Sorathia et al. (2010) were in a "low" 
state can be found in the data of Figure 16. This data of this fig- 
ure suggest that approximately 80% of the total accretion stress 
is found in modes on scales with kyH/2n < 0.3, that is in modes 
with toroidal length scales greater than ~ 3H, corresponding to 
an approximate angular extent of ~ 40°. The simulations de- 
scribed by Fromang & Nelson (2006); Sorathia et al. (2010) 
utilized toroidal domains with angular extents 45° and 60° re- 
spectively. Our results suggest that in the turbulent disk 80% of 
the accretion stress is found on scales greater than 40° and as 
such the simulations described by Fromang & Nelson (2006); 
Sorathia et al. (2010) must significantly underestimate the to- 
tal accretion stress. The observation that the simulations pre- 
sented by these authors exhibit similar vertical stress profiles 
to the "low" state described by Simon et al. (2011) emphasizes 
the non-linear nature of accretion disk turbulence. Choices re- 
garding toroidal domain size can have a profound impact on 
the properties of the turbulence itself. An interesting question is 
therefore whether extending the toroidal domain from n/2 to 
2n will again influence the turbulence. Previous studies suggest 
that this is not the case (Hawley et al. 2001), however, further 
investigation of this issue is necessary. These calculations are in 
progress and will be presented in future work. 



6 STRUCTURE OF THE TURBULENCE 

The results of the preceding section, taken in combination of 
those of Fromang & Nelson (2006); Sorathia et al. (2010), sug- 
gest that the extent of the toroidal domain in global simulations 
plays a key role in determining both the saturation level and the 
vertical profile of accretion stress within these simulations. This 
result is somewhat surprising as the toroidal domain size utilized 
by these authors (7r/4 and n/3 respectively) is sufficient to cap- 
ture many disk scale heights, which should be sufficient since 
results from shearing box calculations suggest that correlations 
in the turbulence are local (Guan et al. 2009). This result points 
to a more fundamental question about the nature of MRI-driven 
MHD turbulence within accretion disks, namely, is the nature 
of the turbulence local (correlations on size scales smaller than 
the disk scale height) or global (correlations on size scales many 
times the disk scale height)? 

Beyond considerations of the influence of azimuthal do- 
main size, there remain questions regarding the behavior of the 
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Figure 18. Two-dimensional power spectrum of magnetic field strength 
(left panel) and kinetic energy (right panel). Each panel shows 
^ I (P^(|/c|)} (solid lines) compared to the power spectrum expected 
for isotropic incompressible turbulence, l/cl"^^/"^ (dashed lines) where 
\k\ = + k^. The power spectra are calculated as described in §3 
and time-averaged over the period AT = 11.5 — 19Porb(^ = l^rg). 
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Figure 19. Power spectra of temporal fluctuations in magnetic field 
strength (left panel) and kinetic energy (right panel). In both panels, 
solid black curves denote results from simulation data, the upper curve 
fluctuations on spatial scales larger than a disk scale height, the lower 
curve fluctuations on spatial scales smaller than a disk scale height. The 
solid red curves represent the result of rebinning the raw power spectra 
onto a grid a factor of 10 coarser in frequency space. The dashed line 
denotes the scaling f~^. The power spectra are calculated as described 
in §3. 



periodic variations in the magnetic field described in §4.1, in- 
cluding: How are these fluctuations arranged spatially? Is there 
an energy injection scale? What temporal variability patterns 
characterize the dynamo? Is there a distinction between global 
and local behavior in this system? In this section, we address the 
question of whether MRI-driven MHD turbulence within accre- 
tion disks is local or global in nature and investigate the prop- 
erties of the turbulence on both small and large scales. Our pri- 
mary diagnostic in performing this analysis are Fourier trans- 
forms of scalar (rather than vector) quantities, calculated as de- 
scribed in §3. We utilize scalar quantities for reasons of both 
simplicity and computational tractability (this approach reduces 
the required number of Fourier transforms by a factor three). 
We have found that while this approximation results in some 
quantitative changes to the outcome of these calculations, the 
qualitative conclusions that we draw remain the same. We note 
that other limitations of our calculations, e.g. finite grid resolu- 
tion and restricted vertical domain will play at least as important 
a role in determining precise quantitative outcomes as our use 
of scalar rather than vector quantities and so, to this extent, we 
argue that this approximation is justified. The remainder of the 
section examines the power spectrum of the magnetic field and 
density fluctuations on the poloidal plane, temporal fluctuations 
of these quantities at small and large spatial scales, fluctuations 
in the magnetic field and density within the disk body and finally 
the autocorrelation functions of these quantities, again within 
the disk body. 



6. 1 Characteristics of Fluctuations in the Magnetic and 
Kinetic Energies 

As a first step in understanding the behavior of MRI-driven MHD 
turbulence within the disk, we further investigate the periodic 
variations in the poloidal structure of the magnetic field dis- 
cussed in §4.1. The first step in this discussion is to under- 



stand how structures in the magnetic and kinetic energies are 
arranged spatially. We probe this aspect of the simulations via 
the shell- and time-averaged two-dimensional Fourier transform 
of toroidally-averaged, spatially normalized simulation data. 



(P|Bp(|fc|)) 



and 



(Pp(ifci)) 



calculated as described 



in §3.2. These data, time-averaged over AT = 10 — 20Porb(^ = 
ISvg) using 3000 frames and calculated from simulation data 
in the region 5 < r/r^ < 15, |Z| < 5 are shown in Figure 18. 
One point of contrast between the two power spectra is that the 
total power in kinetic energy fluctuations is approximately an or- 
der of magnitude less than that in magnetic energy fluctuations. 
Apart from this, the power spectra are remarkably similar. The 
scaling of both power spectra with \k\ is well described by the 
power law consistent with the Kolmogorov spectrum for 

homogeneous incompressible turbulence (Hawley et al. 1995). 
In both cases, the power spectrum is rather featureless, with 
no evidence of a break at large scales, a behavior reminiscent 
of that reported by Hawley et al. (1995); Simon et al. (2011) 
for unstratified shearing boxes. These results suggest that en- 
ergy is injected at the largest scales within the turbulence and 
is transferred to progressively smaller scales via a direct cas- 
cade, a result consistent with those of Lesur & Longaretti (2010). 
That the power spectrum of the magnetic and kinetic energies is 
approximately incompressible at small scales is perhaps unsur- 
prising, as velocity fluctuations and magnetic field strengths are 
largely subthermal in these simulations close to the disk mid- 
plane, |Z| < H. The extension of this same power law to large 
scales therefore implies that turbulent fluctuations within the 
disk remain incompressible at large scales, despite the increase 
in magnetic field strength and turbulent velocity fluctuations as 
one moves to large scales within the disk (see the discussion of 
§5.1 and §5.2). 

One of the underlying assumptions of the Kolmogorov 
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Figure 20. Structure of ^Cp(Ax)^ (left panel), ^C|g|2(Ax)^ (center panel) and ^Cp|5y|2(Ax)y (right panel) where Ax = (Ax, Ay, Az) along the plane 
Az = 0. The correlation functions are calculated over the region 5 < r/r^ < 11, |Z| < ±H) and < <nl2, time- averaged over AT = 16 — 19Porb(^ = 
ISrg) and normalized to their value at Ax = 0. 



model for homogeneous, incompressible turbulence is that the 
fluctuations are self-similar; that is, fluctuations on different 
scales are statistically indistinguishable. One way to probe 
whether this is the case is to examine the character of tempo- 
ral fluctuations in the magnetic and kinetic energies on large 
{\k\H/2n < 1) versus small {\k\H/2n > 1) spatial scales. 
This is accomplished by taking the three-dimensional Fourier 
transform of toroidally-averaged, spatially normalized simula- 



tion data. 



and 



%m,n) 



as described in 



§3.2. We perform a one-dimensional average over \k\H/2n < 1 
and |/c|H/27r > 1 to yield the temporal power spectrum on large 
and small spatial scales respectively (see §3.2). In performing 
these calculations, we utilize simulation data from the same do- 
main as above, 5 < r/r^ < 15, |Z| < 5 and include contributions 
from frames in the range AT = 10 - 20Porb(^ = ISrg). These 
data are shown in Figure 19. Whilst both sets of power spec- 
tra are dominated by contributions from large scales (a conse- 
quence of the |/c|~^^/^ scaling described above), there are now 
clear contrasts between the temporal power spectra for these 
two measures of the turbulence. Firstly, the magnetic energies 
are characterized by broadband power extending over the entire 
frequency range and characterized by the power law, oc f~^^^. 
Kinetic energies are instead characterized by a broken power 
law, cx: at frequencies higher than f P^j^jj^r = 10rs)/2n > 10 
and oc f~^-^ at frequencies lower than this. In addition, we find 
evidence for a break at / ~ Porbi^ = lO^s) the power spectra 
for the magnetic energy at large spatial scales, a feature that 
is absent for the magnetic energy at small spatial scales and 
also in the kinetic energy at all scales. These results suggest 
that while the temporal fluctuations in the kinetic energy are 
similar at different scales within the turbulence, the same is not 
true for the magnetic energy. That is, the temporal behaviour 
of the magnetic field at large versus small scales is not consis- 
tent with the underlying Kolmogorov model for homogeneous, 
incompressible turbulence, namely that fluctuations on different 
scales are statistically indistinguishable. In the next section, we 



examine properties of the temporal correlation function in these 
quantities to confirm this suggestion. 



6.2 Correlation Functions 

Guan et al. (2009) use unstratified local shearing box calcula- 
tions to demonstrate that the magnetic autocorrelation function 
is localized on the x — y (r — 0) plane with a correlation length, 
A, along the major axis of the correlation function of A ~ 0.3H 
for net toroidal magnetic field geometries, suggesting that the 
turbulence is local. The autocorrelation function is useful in ad- 
dressing these questions as it provides an improved statistical 
measure of properties of the turbulence at large spatial scales 
(Guan et al. 2009). Davis et al. (2010) examine the shape of the 
magnetic autocorrelation function on the x — y plane for large 
radial domain, vertically stratified shearing boxes. These authors 
obtain results approximately consistent with those of Guan et al. 
(2009) for zero net flux configurations at small scales, but also 
find the existence of large scale correlations which enforce uni- 
formity in the Maxwell stress and magnetic energies over many 
scale heights, suggesting that global correlations exist within the 
turbulence. That global correlations do exist within the turbu- 
lence was further demonstrated by Nelson & Gressel (2010). 
These authors compare the amplitude of density fluctuations ob- 
tained in vertically stratified shearing boxes with those obtained 
in global simulations, finding that shearing boxes of dimensions 
4H X 16H X 2H in x, y,z are necessary to provide correlations in 
the density which are consistent with those obtained from global 
simulations performed using a 7r/2 toroidal domain. This is be- 
cause at least six scale heights are necessary in the azimuthal do- 
main to correctly capture the excitation of spiral density waves 
Heinemann & Papaloizou (2009a,b). As suggested in §5.1, these 
spiral density waves are intimately tied to the strength of ve- 
locity fluctuations within the disk and hence to the turbulence 
itself. 

If an azimuthal domain size of sixteen scale heights is nec- 
essary to accurately reproduce the results oi n/2 toroidal do- 
main global simulations within shearing boxes, it is worth ask- 
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Figure 21. Slices through the major (solid lines), minor (dashed lines) and vertical (dot-dash lines) axes of (Cp^ (left panel), (Cj^p^ (center panel) 
and ^Cp|5y|2 ^ (right panel), calculated as described in Fig. 20 





Figure 22. Lifetimes of modes with Ax > H (solid lines) compared to modes with Ax < H (dashed lines) for ^ ^Cp(Ax, Aj, At)^ ^ (left panel), 
^ (Cjgp (Ax, Aj, Af )^ ^ (center panel) and ^ ^Cpj^yp (Ax, Ay, At)^ ^ . The correlation functions are calculated from vertically integrated data over the 
region 5 < r/rg < 11, <(/) < n/2 where the vertical integral is performed over |Z| < ±H and over the time-interval At = ±2Porb(r = Sr^). The 
resulting correlation functions are then time-averaged over AT = 16 — 19Porb(^ = l^rg) and normalized to their value at Ax, Ay, At = 0. 



ing whether n/2 toroidal domain global simulations accurately 
reproduce the results of full 2n toroidal domain global simu- 
lations? Hawley et al. (2001) performs an explicit comparison 
of cylindrical global simulations simulations computed using a 
full 2n domain in the toroidal dimension versus those using a 
restricted n/2 domain, finding that, whilst the linear growth 
stage of the toroidal field MRI can be influenced by use of a 
restricted domain, properties of the quasi-steady state are ap- 
proximately similar between the two domains. In particular, the 
one-dimensional toroidal power spectra of the density and mag- 
netic field components were found to be approximately inde- 
pendent of the increase in toroidal domain size from 7r/2 to 271. 
Taken together, these results suggest that it is necessary to con- 
sider toroidal domains at least a factor 16 greater than the disk 
scale height in order to obtain accretion stresses independent of 
the domain size. 

To further investigate these ideas, we have calculated three- 
point auto-correlation functions of the gas density, ^Cp(Ax)^, 

the magnetic energy density, ^C|b|2(Ax)^ and the kinetic energy 
density, ^Cpi^vpCAx)^ (where Ax = (Ax,Ay,Az)) with the 
goal of determining correlation lengths within the turbulence. 
Each of these auto-correlation functions are calculated over the 
region 5 < r/vg < 11 inside the disk body (i.e. |Z| < ±H) 
and over the entire toroidal domain (i.e. < (j) < n/2). The 
extent of the radial domain is chosen so that the correlation 
functions pass through zero before interacting with the radial 
boundaries. The extent of the domain in the vertical dimen- 



Table 1. Correlation Lengths 



Autocorrelation function 










3 


0.38 


1 




2 


0.25 


0.13 




1 


0.25 


0.13 



sion is chosen both for the reasons described in Davis et al. 
(2010), namely so that large scale correlations from the regions 
close to the vertical boundaries do not influence the shape of 
the correlation function and so that we are confident that the 
included region corresponds to incompressible turbulence. The 
correlation functions are calculated as described in §3.2 over 
AT = 16 - 19Porb(^ = ISrg) (corresponding to ~ 10 orbits at 
Sr^) from 204 full three-dimensional data dumps (a resolution 
of ~ 20 dumps per orbital period at Sr^). 

Figure 20 shows each correlation function on the x — y 
plane at zero offset from the vertical axis, AZ = 0. Before dis- 
cussing the properties of these functions in detail, we recall that 
the ratio of the Maxwell stress to the magnetic pressure defines 
a characteristic tilt angle of the magnetic field with respect to 
the toroidal direction (Guan et al. 2009) : 



l(«))l = - 



|B|2 



= 2 sin 9^ cos 0^ = sin 20^ 



(31) 



where 6^ is the 'tilt' angle with respect to the 0-axis. For the sim- 
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ulation presented here, I I ~ 0-3 (see Figure 5) within 

the disk body (|Z| < H), where we are confident that the tur- 
bulence is incompressible, implying that 6^ ~ 9°, somewhat 
smaller than 6^ = 12° reported by Guan et al. (2009) for net 
toroidal field unstratified shearing boxes computed at 32 zones 
per scale height. The data of Figure 20 show that both the auto- 
correlation function in the density and the magnetic field are 
both characterized by ellipses, with the major axis tilted with re- 
spect to the y-axis at an angle consistent with 6^ ~ 9°. The auto- 
correlation function of the kinetic energy density is more circular 
than these two previous measures and as such it is harder to as- 
sess the alignment of the major axis. What evidence there is does 
suggest that the tilt angle of this auto-correlation function is also 
consistent with 6^ ~ 9°. To determine the correlation lengths 
along the major, A^^^, minor, A^^j^^ and vertical axes, A^g^t ? ro- 
tate each correlation function through 6^ ~ 9° on the x—y plane 
and measure the distance from the axis where the correlation 
function falls to of its maximum value. Assuming that the 
correlation functions along each axis take an exponential pro- 
file, (C(Ax)) = (C(0)) e"^^/^ (as in Guan et al. 2009), we can 
then determine the correlation lengths, A along each axis. Plots 
showing the shape of the correlation function along the major, 
minor and vertical axes are shown in Figure 21 and the corre- 
lation lengths are given in Table 1. As these data make clear, 
there is a hierarchy of scales within the turbulence. The longest 
correlation lengths exist in the density, correlations in the mag- 
netic field exist inside these and finally the kinetic energy has the 
smallest correlations lengths. All of these correlation lengths are 
significantly longer (in units of the disk scale height) than those 
reported by Guan et al. (2009) for a net toroidal field model. 
We note, however, that the correlation length along the major 
axis in the density is identical to that found by Nelson & Gressel 
(2010), suggesting that these extended correlation lengths are 
due to the excitement of spiral density waves in the disk (Heine- 
mann & Papaloizou 2009a,b). Finally, the correlation length in 
the density coincides exactly with the break in the toroidal field 
power spectrum found in both the magnetic field and the accre- 
tion stress in §5.2 and §5.3, suggesting an intimate link between 
the the formation of spiral density waves and fluctuations in the 
magnetic field and accretion stresses. We leave detailed investi- 
gation of this possibility to future work. 

The discussion of §6.1 suggests that fluctuations within the 
turbulence can display different temporal variability patterns at 
different spatial scales. If correct, this would indicate that the 
turbulence is not self-similar at different spatial scales. We can 
further investigate this suggestion by examining the lifetimes of 
modes at large versus small spatial scales by using the scale fil- 
tered space-time auto-correlation functions described in §3.2. 
Figure 22 shows the lifetimes of modes in the density, magnetic 
field and perturbed kinetic and large and small spatial scales. For 
each of these quantities, the lifetimes of modes on small spatial 
scales T ~ O.OSPorb (estimated using (C(At)) = (C(0)) e"^'/^) 
where t is the mode lifetime. The lifetimes of modes on large 
spatial scales display different behavior. Fluctuations in the den- 
sity on large spatial scales do not de-correlate on the timescales 
considered here, indicating that they are long-lived compared to 
the local orbital time. Modes on large spatial scales in the mag- 
netic energy density have t ~ O.SPo^bJ whilst those in the kinetic 
energy density have t ^ 0.25Porb- ^e also note that the form 
of the temporal correlation function in the magnetic energy on 
large scales has a different dependence on At at At > O.SP^rb 
to At < O.SPor^; this is the signature of the break in the tempo- 



ral power spectrum in magnetic energy on large spatial scales at 
/ Porb discussed in §6.1. These data show that the lifetimes 
of modes on large scales are significantly longer than those on 
small scales, with typical lifetimes of ~ P^rb in the former case, 
compared to lifetimes ~ OAP^rb in the latter, confirming the sug- 
gestion of §6.1 that the properties of the turbulence are not sta- 
tistically indistinguishable at different spatial scales. 



7 A LOCAL FLUX-STRESS RELATIONSHIP 

In unstratified shearing box simulations, there is a direct re- 
lationship between the strength of the vertical magnetic flux 
threading the domain and the time-averaged accretion stress 
arising from the MRI (Hawley et al. 1995). In a box of fixed 
physical size, the relation takes the simple form (Pessah et al. 
2007), 

-^^-B.^-^MRr (32) 

where A^^^ is the wavelength of the fastest unstable mode of 
the MRI in the vertical direction for the specified net initial mag- 
netic field B^. The proportionality is limited in the high field limit 
by the requirement that the most unstable wavelength must fit 
within the box, if instead A^^^ > L, then the MRI freezes out 
and the accretion stress drops to zero. This is a physical limit 
that would suppress (at least) the linear growth of the MRI in 
a real disk threaded by a field with A^^^ » H. The proportion- 
ality is also limited, purely numerically, in the low field limit by 
the requirement that A^^^ must be resolved, i.e. that A^^^ > A^, 
where A^ is the vertical grid spacing. Because unstratified zero- 
net flux local simulations fail to converge in the absence of phys- 
ical dissipation (Fromang & Papaloizou 2007), a consequence 
for those simulations only is that the mean accretion stress is en- 
tirely determined by the boundary conditions for the magnetic 
field and the grid resolution (Pessah et al. 2007). 

Although the existence of a flux-stress relation has a 
straightforward anchor in the linear physics of the MRI, its rel- 
evance to systems more realistic than the unstratified shearing 
box remains unclear. We first note that stratified shearing boxes 
do not exhibit the same pathological convergence properties as 
unstratified simulations, even when only numerical dissipation 
is present. Rather, such simulations converge to a non-zero ac- 
cretion stress even in the limit of zero net flux (Davis et al. 
2010). This suggests that for local, stratified disks, there ought 
to exist a threshold net flux, below which the properties of the 
disk turbulence are universal. A flux-stress relation would then 
apply only for larger net magnetic fields, up to the (physical) 
limit where MRI quenching occurs. The simulations necessary to 
test such expectations in detail have not yet been done, though 
there is considerable evidence that vertical flux continues to 
stimulate MRI-driven turbulence irrespective of the presence or 
absence of stratification. 

In local simulations, the total vertical flux threading the do- 
main is set by the initial conditions^. Such simulations cannot 
address the question of what would be the self-consistent distri- 
bution of magnetic flux threading different regions of a disk. In 
a real disk, the action of the MRI, together with dynamical pro- 
cesses in the disk corona, may combine to generate local regions 

^ Of course, in a large shearing box the distribution of flux threading the 
mid-plane is a dynamically evolved quantity, and is of physical interest. 
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Figure 23. Time- averaged local correlation between the vertical magnetic flux and the accretion stress within the disk body (|Z| < H) for three different 
radial regions of the disk. The vertical flux is measured in terms of the number of vertical grid cells per fastest unstable mode of the vertical field MRI, 
A]vjK//A and the accretion stress is measured in units of the gas pressure, W^/Pg. The color contours show the fractional distribution of ring volume 
over the parameter space, the solid line the average accretion stress at each ?imri I ^ the dashed line the cumulative fractional volume of the ring at 
each X^^il A (where the numerical value can be read directly from the j-axis). From left-to-right, the panels show data for 3 < r/r^ < 5, 5 < r/r^ < 7 
and 7 <r Irg < 9. The time -averaging is performed over AT = 11.5 — 19Porb(^ = l^rg) from 3000 complete three-dimensional data dumps. 



of strong flux, even if the disk as a whole has vanishing net field. 
Conversely, vertical flux in a disk initially threaded by a net field 
way well be transported over large radial distances, or expelled 
from the disk entirely (Beckwith et al. 2009). Global simulations 
are needed to address these possibilities. Two questions are of 
particular import. First, what is the form of (and dispersion in) 
the flux-stress relation measured globally, and on what spatial 
scales does it apply? It is necessary to check whether the rela- 
tion is causal - in the expected sense of the flux influencing the 
stress rather than the other way around - since this is not nec- 
essarily true globally. Second, is the distribution of vertical flux 
threading the disk sufficient as to boost the accretion stress, over 
and beyond the expectations of a model in which the disk has 
zero net flux on all relevant local scales? It is possible to imag- 
ine a situation in which the self-consistent coronal field has a 
dominant feedback effect on the dynamics at the disk mid-plane 
- increasing the accretion stress from the ((a)) ~ 10~^ values 
measured in many simulations to values closer to those observa- 
tionally inferred - but the existence of such a regime has yet to 
be demonstrated. 

Sorathia et al. (2010) made an initial examination of these 
issues, using data from a series of global disk simulations that 
were initialized with small loops of poloidal field. They analyzed 
the model disk by considering the flux and stress present in co- 
moving patches of a single, fixed spatial scale, concluding both 
that a causal flux-stress relation existed and that the distribution 
of vertical flux was strong enough as to be dynamically interest- 
ing. Here, we revisit these questions. We take advantage of the 
fact that our simulation is initialized with a net toroidal field, 
which means that any flux-stress relationship within this simu- 
lation must arise due to vertical magnetic field generated from 
turbulence arising from the MRI. Compared to the simulations 
analyzed by Sorathia et al. (2010), we also attain substantially 
better resolution (in the final, saturated state) of the most un- 
stable linear MRI modes. This permits a separation between nu- 



merical effects that are known to occur on the grid scale, and 
the physical effects at larger scales that are of primary interest. 

To assess whether or not a flux-stress relationship exists 
in our simulation, we first determine the fractional volume of 
the disk that is instantaneously threaded by a given vertical flux 
and that has a given accretion stress. By working with this two- 
dimensional distribution as our fundamental quantity, we avoid 
the need to make an essentially arbitrary choice of spatial scale 
over which to average the flux and the stress. Operationally, 
we divide the disk into three rings spanning 5 < r/r^ < 7, 
7 < r /rg < 9 and 9 < r/r^ < 11. Next, we utilize calculate 
the ratio of the wavelength of the fastest unstable mode of the 
vertical field MRI, A^^^ to the physical grid spacing in the ver- 
tical direction, for each cell within the disk body (|Z| < H), 
which we denote by = A^^^A^. We utilize this same data 
to calculate the magnetic (Maxwell) accretion stress in units of 
the gas pressure, = —W^^jP^, again for each cell within the 
disk body. We use a volume weighted binning procedure to cre- 
ate a distribution function for each ring describing the relation- 
ship between vertical flux and accretion stress, which we nor- 
malize to the total volume of the ring. This procedure yields a 
volume weighted distribution function describing the relation- 
ship between vertical flux and accretion stress, which we denote 
by/(F^a-): 



rlTl ^nF\ an dF^ dor 

J j j dVdF^ da^ 



(33) 



where 5V(F^,a^) is the volume element of a cell threaded by 

vertical flux F^ and an accretion stress and j j j dVdF^da^ 
is the total ring volume multiplied by the area of the (F^,a^) 
plane. We can use this distribution function to find the mean 
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Figure 24. Integrated cross-correlation, ^ ^ ^C^z ^v^m (At)y y y calcu- 
lated as described in §3.2 time- averaged between orbits 10 — 12 (top 
left), 12 - 14.5 (top right), 14.5 - 16 (bottom left) and 16 - 19 (bot- 
tom right) where the orbital period is measured at r = 15r5. The 

cross-correlation function is defined such that if ^ ^ {^^f^ ^'^'^'^) ) ) 
is maximized at negative (positive) At, then fluctuations in F^{r, (p, t) 
lead (trail) fluctuations in W^(r, cj), t). 



accretion stress, associated with a given vertical flux: 



{5V(Fn} : 



action of the i 
:ical flux: 



(34) 



and also the fraction of the ring volume, {5V} that is threaded 
by a given vertical flux: 



jdc 



(35) 



These data are shown in Figure 23, time-averaged over AT = 
11.5-19Porb(^ = ISr^) using 20 dumps per ISCO orbit (approx- 
imately 3000 dumps in all). In these figures, the color contours 
show the volume-weighted distribution function, /(F^, a^), the 
dashed white line the fraction of the ring volume threaded by 
a given vertical flux, {5V(F^)} (where {5V} can be read from 
the y-axis scale) and the solid white line the mean accretion 
accretion stress associated with a given vertical flux, |a^(F^) j. 

The data of Figure 23 demonstrates that each of the rings 
listed above exhibit flux-stress relationships that are broadly 
similar. The majority of each ring is threaded by significant ver- 
tical flux, > 1, a result we have confirmed by direct inspec- 
tion of simulation data. There is a broad range (3 — 4 orders 
of magnitude) of accretion stresses associated with a given F^. 
Nevertheless, the mean accretion stress associated with a given 
vertical flux, |a^(F^) j displays a behavior that is broadly con- 
sistent with the expectation described above, namely that for 
weak vertical fields threading a cell, the mean accretion stress is 
approximately independent of the strength of the vertical field, 
whilst for strong vertical fields threading a cell, the mean accre- 
tion stress is roughly proportional to the strength of the vertical 
field. Notably, the transition between these two regimes takes 
place where the wavelength of the vertical field MRI is resolved 
by approximately 5 cells, consistent with the results of Hawley 
et al. (1995), where it was found that approximately this num- 
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Figure 25. Structure of {^{^Cpz yY^^(Ax,Ay,At)j j at At = calcu- 
lated as described in §3.2 and time- averaged between orbits 10 — 12 
(top left), 12 - 14.5 (top right), 14.5 - 16 (bottom left) and 16 - 19 
(bottom right) where the orbital period is measured at r = 15r5. 



ber of cells per fastest unstable mode was necessary to repro- 
duce numerically the expectation for the linear growth rate aris- 
ing from analytic theory. This result in particular gives us confi- 
dence that there is a physical relationship between vertical flux 
and accretion stress, as we do not need to resort to arguments 
regarding the growth rate of long wavelength modes (Sorathia 
et al. 2010). 

The preceding discussion suggests that there is a local re- 
lationship between vertical flux and accretion stress operating 
within the disk. It does not, however, describe the sense of that 
relationship, i.e does a fluctuation in the vertical flux lead to 
a fluctuation in the accretion stress or vice versa? Determining 
which of these possibilities is the case is clearly crucial in de- 
ciding whether or not this flux-stress relationship is physical. 
Sorathia et al. (2010) use co-moving wedges to calculate the 
temporal correlation between vertical flux and accretion stress, 
where the azimuthal velocity of the fluid is used to track a given 
co-moving wedge between timesteps. This approach however, 
assumes that fluctuations in the magnetic field propagate purely 
in the azimuthal direction with the fluid rotation velocity, which 
may not be appropriate in a turbulent magnetized disk. We in- 
stead address the causal relationship between vertical flux and 
accretion stress using the three-point cross-correlation function 
approach outlined in §3.2, an approach which bypasses the need 
to consider a co-moving wedge as it calculates the correlation 
between all points within the domain of the correlation func- 
tion simultaneously. The cross-correlation functions are calcu- 
lated over the radial region S < r /r^ < 11 using vertically inte- 
grated data within the accretion disk body (|Z| < H) for 2 orbits 
at the center of the radial domain, r = Sr^ using 250 frames. 
We time-average the cross-correlation functions at four different 
times within the evolution, chosen to coincide with oscillations 
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in the dynamo cycle discussed in §4.1, namely orbits 10 — 12, 
12 - 14.5 14.5 - 16 and 16-19 where the orbital period is 
measured at r = 15r5. Each time-average is computed using 70 
dumps per orbit at r = 15r5. 

Figure 24 shows the total amplitude of the cross-correlation 
function as a function of the orbital time at Sr^ for each of these 
different periods in the dynamo cycle. Here, biasing of the corre- 
lation function to negative (positive) At/P^rb indicates that fluc- 
tuations in the vertical flux (accretion stress) lead fluctuations 
in accretion stress (vertical flux). Clearly, at orbits 12 — 14.5 and 
14.5 — 16, we see that fluctuations in the vertical flux lead fluctu- 
ations in the accretion stress by approximately 0.2Porb(^ = Sr^), 
a period consistent with the results of Simon et al. (2009) for 
the typical lifetime of a turbulent fluctuation. The causal nature 
of the relationship between vertical flux and accretion stress at 
orbits 10 — 12 and 16 — 19 is less clear. Here, the correlation 
functions are double peaked, with the peaks located on either 
side of At = 0. This result suggest that at these points in the dy- 
namo cycle, fluctuations in the vertical flux lead to fluctuations 
in the toroidal (or radial) magnetic field, which results in a fluc- 
tuation in the accretion stress. The toroidal (or radial) magnetic 
field fluctuation then leads to a new vertical flux fluctuation 
some time At ~ 0.6Porb(^ = Sr^) later, a process reminiscent 
of the simple model for a magnetic dynamo in an accretion disk 
described by Tout & Pringle (1992). These results would sug- 
gest then that during orbits 10 — 12 and 16 — 19, we should see 
rapid rearrangement of the magnetic field as vertical, toroidal 
(and presumably radial fields) couple together, whilst during 
orbits 12 — 14.5 and 14.5 — 16, the structure of the magnetic 
field should be more stable. Figure 4 suggests that this is indeed 
the case, at orbits 12 — 14.5 and 14.5 — 16, we see a relatively 
ordered toroidal magnetic field within the simulation domain, 
whilst at orbits 10 — 12 and 16—19, the toroidal magnetic field 
is somewhat less structured. 

A final probe of the flux-stress relationship within the 
disk body is the shape of the cross-correlation function on the 
(Ax, Ay)-plane at zero temporal offset. At = 0, calculated from 
the same data as used to create the temporal correlation func- 
tions shown in Figure 24. These data are shown in Figure 25, 
again for 10-12, 12-14.5 14.5 -16 and 16-19 where the or- 
bital period is measured at r = 15r5. Comparing the data of this 
figure with that of Figure 20, we see that the flux-stress cross- 
correlation function has approximately the same tilt angle with 
respect to the Ay axis as found for ^Cp(Ax)^, etc and that the 
correlation length along the major axis lies approximately inside 
that for fluctuations in the density. Whilst the autocorrelation 
function in a quantity must be symmetric with its maxima lo- 
cated at zero-offset by definition, no such requirement exists for 
the cross-correlation function between two quantities. We find 
that the cross-correlation function between vertical flux and ac- 
cretion stress is double peaked on the (Ax, Ay)-plane, with the 
peaks lying along the major axis of the correlation function cen- 
tered on |Ax| 0.25H. At orbits 10 - 12 and 14.5 - 16, the ra- 
dially outer peak is greatest in amplitude, whilst at 14.5 — 16 the 
peaks are approximately equal in amplitude and orbits 16—19, 
the radially inner peak is greatest in amplitude. Overall, this sug- 
gests that there are vertical fieldlines penetrating the disk body 
which link together adjacent radii in a manner reminiscent of 
that suggested by Tout & Pringle (1996). Combining this result 
with those discussed in the preceding paragraph suggests that as 
toroidal field rises out of the disk body (presumably due to mag- 
netic buoyancy), vertical field is created that penetrates the disk 



midplane. If this vertical field is well resolved in terms of the 
fastest unstable mode of the vertical field MRI, then we measure 
an enhanced accretion stress associated with the operation of 
the vertical field MRI which also creates new toroidal and radial 
magnetic fields within the disk. These new magnetic fields even- 
tually become buoyantly unstable and the process repeats. This 
process is strongly reminiscent of the dynamo model described 
by Tout & Pringle (1992), which one would hope could provide 
a detailed analytic framework in which to describe the results 
presented here. We leave such calculations to future work. 



8 SUMMARY, DISCUSSION AND CONCLUSIONS 

Global simulations of magnetized thin accretion disks are re- 
quired to study the physics of the MRI on large scales (both spa- 
tial and temporal), and to make predictions for the structure of 
turbulent accretion disks that can be tested observationally. To 
be useful, such calculations must ideally be computed at resolu- 
tions which approach that of local simulations, which has only 
recently become feasible. Here, we have presented an analysis of 
a simulation of a global, magnetized, thin (H/R ~ 0.07) accre- 
tion disk designed to investigate the properties of MRI-driven 
turbulence. The simulation was initialized with a moderately 
strong net toroidal field (similar to, e.g. Hawley & Krolik 2002; 
Fromang & Nelson 2006; Beckwith et al. 2008a), but it rapidly 
loses memory of its initial conditions and reverts to a state that is 
consistent with models initialized with zero net flux. The com- 
putation used a second order Godunov scheme with accurate 
fluxes at a poloidal resolution comparable to moderately well- 
resolved local simulations (see e.g. Simon et al. 2011). Our al- 
gorithmic choices have been shown to capture the linear growth 
stage of the MRI accurately (Flock et al. 2010), and likely yield 
improved accuracy at fixed spatial resolution over prior simula- 
tions. The results allow us to make a quantitative assessment of 
the structure and locality of the resulting turbulence, and inform 
a qualitative discussion of the implications for observations and 
simplified models of disk dynamos. 

Our results for the locality of MRI-driven disk turbulence 
suggest a nuanced picture, in which some aspects of the turbu- 
lence are well-described by a local model, whereas others re- 
quire a global treatment. From an analysis of the spatial two- 
point correlation functions, we find that accretion disk turbu- 
lence, whilst subsonic, contains significant correlations on scales 
> H. The longest correlation lengths exist within the density, for 
which Xj^aj = 3H, followed by the magnetic energy iX^^aj — 2H) 
and then the kinetic energy (X^^aj = ^H), implying that the 
largest scales within the turbulence are controlled by the density. 
Prominent spiral density waves are observed in these simula- 
tions and the correlation length along the major axis of the den- 
sity correlation function suggests that it is these structures that 
set the size of turbulent fluctuations within the disk (see e.g. Nel- 
son & Gressel 2010). If so, then this implies that correctly captur- 
ing the formation of these structures is essential to understand 
the properties of turbulence within the disk. Across a range of 
spatial scales, the fluctuations in the magnetic and kinetic en- 
ergies on the poloidal plane are arranged in a fashion consis- 
tent with expectations arising from homogeneous isotropic tur- 
bulence (i.e. see e.g. Hawley et al. 1995). 

Although the spatial power spectrum of the turbulent fields 
is consistent with a simple incompressible turbulence model, the 
temporal behavior evidences greater complexity. At the most ba- 
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sic level, structures on large scales within the turbulence have 
lifetimes significantly longer than structures on small scales, as 
one might expect. However, we also observe that the structure 
of the variability varies significantly with spatial scale, break- 
ing the self-similar assumption that underlies simple turbulence 
models. We find that temporal fluctuations in the magnetic field 
exhibit different properties at large Qk\H/2n < 1) versus small 
(\k\H/2n > 1) spatial scales. This is in contrast to the ki- 
netic energy, where similar temporal power spectra are found 
at both spatial scales. We tentatively attribute this behavior to 
the lack of a clean separation between the energy injection scale 
and the dissipation scale, which can lead to a non-local (in k 
space) transfer of energy in MHD turbulence (Lesur & Longaretti 
2010). Unfortunately, determining robustly the range of scales 
over which the MRI taps the shear energy of the disk is not pos- 
sible given our resolution, or with any resolution feasibly attain- 
able in a global calculation. Large shearing boxes remain the 
best numerical setups for studying such questions. 

For comparison with real systems, the most basic diagnostic 
of the properties of accretion disk turbulence is the magnitude 
of the time-averaged accretion stress, which can be inferred ob- 
servationally from modeling of thermally unstable disks in dwarf 
novae and X-ray binaries. For systems where the disk is gas pres- 
sure dominated, and hot enough that ideal MHD applies, obser- 
vational estimates suggest ((a)) ^0.1 (Lasota 2001; King et al. 
2007). We measure a value ((a)) ^ 2.5 x 10~^ from our sim- 
ulations, that is larger than that derived from most prior cal- 
culations without net vertical field, but still formally inconsis- 
tent with observations. In principle, this discrepancy could point 
to a physical effect (disks in binaries could be threaded by, or 
spontaneously develop, net vertical fields), but it could also be 
a numerical artifact (the convergence of global simulations has 
not been demonstrated), or be due to a flawed comparison be- 
tween simulations and models of outbursting disks computed 
using classical disk theory. Further work is needed to address 
each of these possibilities. The vertical distribution of the ac- 
cretion stress is approximately constant within ~ 2H of the disk 
midplane. This is consistent with the results of well-resolved ver- 
tically stratified shearing box models (Simon et al. 2011), but is 
in contrast to previous global simulations (see e.g. Fromang & 
Nelson 2006; Sorathia et al. 2010) where a double peaked stress 
profile was measured. We have found that approximately 80% 
of the total accretion stress is located at toroidal angular scales 
> 40°. We attribute the contrasts in the total accretion stress 
and vertical stress distribution between this work and Fromang 
& Nelson (2006); Sorathia et al. (2010) as being due the use 
of toroidal angular domains of 7r/4 and 7r/3 by these authors 
respectively. This leads us to regard the use of toroidal domains 
of extent n/2 essential in order to correctly capture the physics 
of angular momentum transport driven by the MRI. 

Measurements of the turbulent velocity field in disks 
(Horne et al. 1994; Carr et al. 2004) are potentially more power- 
ful probes of disk physics than single-point comparisons of mea- 
sured and simulated accretion stress, provided that a separa- 
tion of turbulent motion from other non-Keplerian flow is possi- 
ble. Protoplanetary disks currently represent the most promising 
observational targets (Hughes et al. 2011). In the ideal MHD 
limit, we find that the vertical distribution of velocity fluctu- 
ations steepens from ~ O.lc^ in the disk midplane to ~ 0.20^ 
in the corona. This implies that magnetized turbulence within 
disks is characterized by turbulent line widths that are between 
0.1 and 0.2 of the local sound speed. The fluctuation amplitude 



that we observe is consistent with previous calculations of MHD 
turbulence in protoplanetary disks Fromang & Nelson (2006). It 
also matches recent observations of such systems Hughes et al. 
(2011), although the importance of non-ideal MHD effects for 
protoplanetary disks means that a quantitative comparison re- 
quires more realistic simulation work than that presented here. 

Largely for reasons of numerical convenience, our simula- 
tion utilized a pseudo-Newtonian potential that results in an in- 
nermost stable circular orbit near the inner boundary of the disk. 
We studied the structure of the disk near and within the ISCO, 
whose detailed properties are important for observational at- 
tempts to measure the spin of accreting black holes (Zhang et al. 
1997; Brenneman & Reynolds 2006; Done et al. 2007; Steiner 
et al. 2010). We found consistency between the measured radial 
turbulent accretion stress distribution within the disk, and the 
expectations of models that assume a stress free condition at the 
ISCO. At and inside the ISCO, accretion stresses are due to large 
scale correlations in radial and toroidal magnetic fields and the 
flow dynamics is controlled by 'flux-freezing' rather than turbu- 
lence. Our results support a picture in which, for thin disks, it is 
the level of net vertical flux at the ISCO that is crucial in deter- 
mining the stress levels there. If the net vertical flux is small or 
zero, then stresses are negligible, whereas significant non-zero 
vertical flux is associated with stresses that could be of obser- 
vational importance (Agol & Krolik 2000). A number of sim- 
ulations - both fully relativistic and pseudo-Newtonian - lend 
credence to this scenario (Reynolds & Armitage 2001; Beckwith 
et al. 2008b; Penna et al. 2010; Noble et al. 2010). 

Numerical simulations - whether they be local or global 
- cannot follow disks over the very long timescales that are 
characteristic of many interesting observational phenomena. It 
is therefore important to understand whether there are features 
of simulated disks that can be abstracted for use in simpler mod- 
els of disk evolution or disk dynamos. One interesting question 
is whether or not a local connection between vertical flux and 
accretion stress (Hawley et al. 1995; Pessah et al. 2007) per- 
sists in global simulations (Sorathia et al. 2010). We find that 
such a relationship does exist within the turbulence and that 
provided the vertical flux threading a given cell is sufficiently 
strong (here determined by the criterion that the wavelength of 
the fastest unstable mode of the vertical field MRI associated 
with the flux threading a given cell is resolved), the vertical flux 
acts as an accurate predictor of the accretion stress. We further 
find, that in a time-average sense, the majority of the disk body 
is threaded by vertical fluxes that are well-resolved by this crite- 
rion. By use of two-point space-time cross-correlation functions 
between the vertical flux and the accretion stress, we find that 
causal sense of the flux-stress relationship depends on which 
point in the dynamo cycle the correlation function is calculated. 
When the toroidal field is well-ordered on the poloidal plane, 
there is a causal connection between vertical flux and accre- 
tion stress. Less ordered toroidal field configurations are associ- 
ated with causal vertical flux-accretion stress that are less well- 
defined. Speculatively, this behavior suggests that as toroidal 
field emerges from the disk body during the dynamo cycle, verti- 
cal fields thread the disk body for sufficient lengths of time that 
they become unstable to the vertical field MRI and thereby de- 
termine the local accretion stress. This process is reminiscent of 
the dynamo model described by Tout & Pringle (1992). 
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